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INLEIDING 


In  dit  rapport  wordt  het  simuleren  van  de  inwendige  ballistische  processen  tijdens  een 
schot  van  een  Otomelara  76/62  behandeld.  Deze  simulaties  zijn  uitge\'oer(l  met  de  variant 
’Grainburning’  van  het  programma  PISCES  2DELK. 

In  het  eerste  deel  zal  kort  worden  ingegaan  op  dit  programma.  Vervolgens  ziillen  de 
gegevens  van  de  Otomelara  worden  behandeld,  waarna  met  enige  resuiiaten  wordt 
afgesloten. 


2  GRAINBURNING  -  PISCES  2DELK 

2.1  Inleiding 

Het  PISCES  2DELK  computer  programma  lost  de  vergelijkingen  op  die  voortkomen  uit  de 
behoudswetten  van  materie,  energie  en  impuls  voor  vloeistoffen  en  vaste  stoffen  in  twee- 
dimensionale  problemen  met  een  axiale  of  axiale  symmetric.  Typische  toepassingen  zijn 
berekeningen  op  het  gebied  van  dynamische  vloeistof  stromingen,  .structuur  dynamica  of 
combinaties  hiervan, 

De  mechanica  van  verschillende  processen  kunnen  in  het  programma  worden  beschreven 
met  een  van  de  vier  volgende  referentie  kaders  (of  een  combinatie  daarvan); 

De  Euler  processor  werkt  met  een  ruimtelijk  constant  niet-deformeerbaar  grid  waarin 
massa,  energie  en  impuls  van  grid-cel  naar  grid-cel  kunnen  bewegen. 

In  de  Lagrange  processor  wordt  het  grid  gedefinieerd  op  het  materiaal.  De  massa  kan  niet 
van  een  grid-cel  in  een  andere  vloeien  en  iiel  grid  kan  uaaruni  vervormen. 

De  thin  shell  processor  is  gebaseerd  op  het  Lagrange  principe,  alleen  in  een  ’dunne 
schaal’  formulering. 

De  rigid  body,  tenslotte,  beschrijft  inerte  lichamen  zonder  interne  dynamica. 

De  Lagrange  en  Euler  formuleringen  zijn  theoretisch  equivalent,  maar  in  de  praktijk 
worden  de  toepassingsgebieden  van  beide  beperkt  door  de  numerieke  benadering, 
waardoor  verschillen  ontstaan.  Grofweg  kan  men  stellen  dat  processen  in  vloeistoffen  en 
gassen  het  best  gesimuleerd  worden  met  de  Euler  processor,  terwijl  voor  deformaties  in 
vaste  stoffen  de  Lagrange  methode  het  geschiktst  is. 

Meerdere  processoren  kunnen  onderling  gekoppeld  worden,  om  bijvoorbeeld  een  systeem 
bestaande  uit  vaste  stoffen  en  vloeistoffen  te  benaderen. 

Omdat  de  te  beschrijven  processen  zo  dynamisch  zijn,  bepaalt  het  programma  zelf  zijn 
tijdstap.  Dit  gebeurt  zodanig  dat  schokgolven  of  trillingen  nooit  in  een  cycle  een  hele 
gridcel  kunnen  passeren.  Merk  op  dat  daardoor  naast  de  golfsnelheid  de  afmeting  van  de 
gridcellen  (met  name  de  kleinste  gridcel)  mede  bepalend  zijn  voor  de  tijdstap. 

Een  speciale  applicatie  van  PISCES  2DELK  is  de  Grain  burning  module,  waarin  het 
mogelijk  is  om  de  overdracht  van  massa,  energie  en  impuls  van  twee  verschillende  (fases 
van)  stoffen  te  beschrijven.  Een  voorbeeld  hiervan  is  het  verbranden  van  kruitkorrels  (en 
de  daarbij  horende  vorming  van  kruitgassen)  in  de  loop  van  een  kanon  tijdens  een  schot. 

Een  uitgebreid  overzicht  van  de  theoretische  structuur  van  PISCES  2DELK  en  de 
applicatie  Grain  burning  is  te  vinden  in  respectievelijk:  de  PISCES  2DELK  theoretical 
manual  (1987)  van  S.  Hancock  [1]  en  het  Technical  Report  TR  87-2118/1  ’Equations  and 
numerical  solution  of  gas-panicle  flow.  GRAINS  level  26’  (1987)  van  P.H.L.  Groenen- 
boom  [2].  Beide  rapporten  komen,  evenals  het  programma  zelf,  van  ’The  MacNeal- 
Schwendler  Company  BV. 
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2.2  In-  en  uitvoer 


De  gegevens  die  het  model  van  het  te  beschrijven  problecm  bevatten,  worden  bet 
programma  op  twee  manieren  aangeboden,  namelijk  via  een  input-file  (extensie  *.I2D)  en 
via  externe  subroutines  die  aan  het  programma  ’gelinked’  worden.  Met  name  dit  laatste 
biedt  de  gebruiker  de  mogelijkheid  om  zijn  model  haast  ongelimiteerd  te  verfijnen  en 
maakt  het  programma  bijzonder  flexibel. 

Er  zijn  verschillende  manieren  van  uitvoer  mogelijk.  In  de  eerste  plaats  kan  een  leesbare 
’print  output-file’  (extensie  *.PRT)  worden  aangemaakt.  Omdat  zo’n  print  output  snel 
enorm  omvangrijk  wordt,  is  het  aan  te  bevelen  om  ook  een  zogenaamde  ’archive-file’ 
(extensie  *.ARC)  te  maken,  die  niet  direct  leesbaar  is  maar  waaruit  met  behulp  van  de 
post-processor  ’DDARCH’  snel  een  selectie  kan  worden  gemaakt  die  wel  leesbaar  is.  Ook 
kunnen  met  behulp  van  de  ’DDPLOT’  en  ’PIPLOT’  post-processoren  uit  de  archive-file 
figuren  (zoals  grafieken,  profielen  en  contour-  en  toestandsdiagrammen)  worden  gemaakt. 
Deze  post-processoren  zijn  eenvoudig  en  snel  te  gebruiken. 

Het  is  ook  mogelijk  om  het  hoofdprogramma  direct  een  figuur  te  laten  maken.  Vaak  is  dit 
echter  moeilijk,  omdat  bijvoorbeeld  v66r  de  berekening  niet  altijd  de  geschikte  tijdstippen 
of  afmetingen  bekend  zijn. 

Het  programma  biedt  de  mogelijkheid  om  een  ’restart-file’  (extensie  *.RES)  te  maken, 
waarin  alle  relevante  gegevens  zoals  die  op  van  te  voren  te  definieren  tijdstippen  (of  cy¬ 
cles)  en/of  aan  het  einde  van  de  berekening  waren  zijn  opgeslagen.  Met  deze  file  kan  men 
de  berekening  laten  doorgaan  op  een  van  deze  momenten. 

Op  het  moment  is  er  (nog)  geen  interactieve  pre-  en  post-processor  bij  het  programma 
beschikbaar.  Dit  is  een  duidelijke  handicap. 


2.3  De  input-file 

In  deze  paragraaf  zal  de  input-file  worden  bekeken.  Omdat  de  manual  die  bij  het  program¬ 
ma  wordt  geleverd  bijzonder  uitgebreid  is,  zal  dit  beknopt  worden  gehouden. 

De  input-file  is  in  twee  stukken  te  verdelen.  Het  eerste  stuk,  in  de  manual  ’global  input’ 
genoemd,  bevat  gegevens  die  van  belang  zijn  voor  de  berekening  in  het  algemeen,  zoals 
bijvoorbeeld  limieten  voor  de  duur  van  het  programma,  de  tijdstapgrootte  en  materiaalge- 
gevens  van  de  stoffen  die  in  het  probleem  voorkomen.  In  het  tweede  stuk,  de  ’local  input’ 
genoemd,  worden  de  subgrids  gedefinieerd,  zowel  wat  betreft  hun  type,  hun  afmeting,  de 
manier  waarop  ze  aan  elkaar  (of  de  ruimte)  zijn  gekoppeld  en  de  wijze  waarop  ze  met 
materiaal  gevuld  zijn.  Zowel  in  het  ’global’  als  in  het  ’local’  gedeelte  zijn  aanwijzingen 
over  de  vorm,  frequentie  en  uitgebreidheid  van  de  output  te  plaatsen. 

Een  voorbeeld  van  een  input-file  is  te  vinden  in  appendix  1,  waar  de  input  voor  het 
Otomelara  76/62  model  is  gegeven.  Hier  is  overigens  te  zien  dat  de  scheiding  in  een 
’global’  en  een  ’local’  gedeelte  op  het  eerste  gezicht  niet  op  te  merken  is,  de  scheiding 
moet  aan  de  kaarten  herkend  worden. 

Het  programma,  althans  de  Grain  burning  variant,  heeft  een  format  gebonden  input,  de 
informatie  dient  op  regels  van  tien  plaahsen  van  acht  characters  gegeven  te  worden. 
Daarbij  moeten  de  gegevens  op  een  vaste,  door  het  programma  bepaalde,  plaats  worden 
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gegeven. 

Waarschijnlijk  heeft  dc  basisvorm  van  PISCES  2DELK  (in  tegenslelling  tot  Grain  binning) 
inmiddels  een  ’free  format’  input.  Dat  dit  echter  nog  niet  het  comfort  van  ecu  interaciieve 
pre-  en  post-processor  met  zich  meebrengt  zal  duidelijk  zijn. 

Enige  speciale  aandacht  verdient  de  manier  waarop  de  materiaaleigenschappen  in  het 
’global’  gedeelte  van  de  input  worden  gegeven.  In  de  MATERIAL  kaart  nioet  de 
toestandsvergelijking  (equation  of  state)  worden  aangeduid.  Het  programma  kent  standaard 
enige  modellen:  de  algemene  gaswet  (gamma  wet),  een  polynoom  vergelijking,  een 
schokgolf  model,  het  Tillotson  model  en  het  Jones-Wilkins-Lee  model.  Daarnaast  is  het 
mogelijk  om  met  behulp  van  een  zelf  te  schrijven  exteme  subroutine  een  eigen  model  te 
gebruiken. 

Enige  regels  verderop  worden  (alleen  voor  de  Grain  burning  variant)  in  de  GRNMISC 
kaarten  enige  materiaalparameters  gegeven  (in  de  COMMENT  kaarten  meteen  daarboven 
wordt  hun  betekenis  gegeven).  Veel  van  deze  gegevens  zijn  echter  niet  met  grote 
nauwkeurigheid  bekend,  hetgeen,  zoals  verderop  uitgebreider  aan  de  orde  zal  komen,  voor 
ons  een  van  de  grote  moeilijkheden  bij  het  simulatie  proces  heeft  gevormd.  Bijvoorbeeld 
over  nummer  vier,  de  proportionality  constant  interfacial  drag  tussen  kruitkorrel  en 
kruitgas,  is  nagenoeg  geen  kwantitatieve  informatie  beschikbaar.  Nu  heeft  de  waarde 
hiervan  gelukkig  geen  dramatische  invloed  op  het  resultaat,  zodat  een  min  of  meer 
acceptabele  en  werkbare  waarde  is  genomen.  Van  veel  groter  belang  zijn  de  burning  rate 
parameters  (nummers  II  en  12),  afkomstig  uif  de  wet  van  Vieille.  Zij  beinvloeden  de 
uitkomsten  zodanig,  dat  zij  tot  op  enkele  procenten  nauwkeurig  bekend  moeten  zijn  om 
een  realistisch  beeld  te  geven.  Dat  dit  problematisch  was  en  dat  waarschijnlijk  een  via  een 
exteme  subroutine  gedefinieerde  altematieve  verbrandingswet  geschikter  is,  wordt  later 
behandeld. 

Tenslotte  definieert  men  in  de  input-file  de  geometrie  van  het  grid  en  de  gewenste 
uitvoerparameters. 


2.4  De  exteme  (user)  subroutines 

Het  programma  biedt  de  mogelijkheid  om  aan  het  basismodel  eigen  varianten  of  aanvullin- 
gen  toe  te  voegen.  Dit  gebeurt  met  behulp  van  zelf  te  schrijven  exteme  subroutines,  die 
aan  het  basisprogramma  ’gelinked’  moeten  worden.  Deze  faciliteit  maakt  het  programma 
bijzonder  plooibaar  en  is  waarschijnlijk  een  van  haar  sterke  punten  in  vergelijking  met 
haar  concurrenten. 

In  appendix  2  is  als  voorbeeld  het  pakket  subroutines  van  het  Otomelara  76/62  probleem 
gegeven.  Het  is  geschreven  in  FORTRAN,  net  zoals  het  programma.  Merk  op  dat  de 
eerste  33  regels  formeel  het  hoofdprogramma  vormen  en  dat  het  werkelijke  standaardpro- 
gramma  als  subroutine  CODE  wordt  aangeroepen  (CALL  CODE).  Waar  in  dit  rapport 
over  hoofd-  of  basisprogramma  wordt  gesproken,  wordt  deze  subroutine  CODE  bedoeld. 

'We  zullen  de  subroutines  nu  nader  bekijken.  Meteen  valt  op  dat  deze  (kunnen)  beginnen 
met  een  uitgebreid  COMMON-blok,  dat  het  mogelijk  maakt  om  variabelen  uit  het 
hoofdprogramma  te  gebruiken. 

De  eerste  subroutine  is  EXFLW2.  Deze  moet  in  de  input-file  (*.12D)  opgeroepen  worden 
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in  de  BOUNDERY  FLOW  kaart  en  bevat  informatie  over  de  wijze  waarop  niateriaal  aan 
de  randen  een  subgrid  kunnen  in-  of  uitstromen.  Hier  is  hij  gebruikt  om  de  aanvuurlading 
in  de  vlampijp  te  simuleren.  De  processen  in  de  vlampijp  zelf  die  besiaat  uit  een  rij 
ongebruikte  gridcellen  worden  niet  gesimuleerd.  Aan  de  randen  van  de  vlampijp  wordt  een 
polygon  gedefinieerd  waaruit  onder  de  aanname  van  kritische  stroining  de  hete  verbran- 
dingsgassen  van  de  ontstekingslading  de  huls  instromen.  De  gaatjes  in  de  vlampijp  worden 
gesimuleerd  middels  een  porositeitsfactor. 

Vervolgens  komt  de  EXEOTW  subroutine,  die  de  toestandsvergelijking  van  de  kruitgassen 
definieert.  In  dit  geval  is  daar  de  covolume  vergelijking  (ook  wel  bekend  als  Noble-Abel 
vergelijking)  voor  gebruikt.  Deze  subroutine  hoeft  niet  gedeclareerd  te  worden  in  de  input- 
file,  maar  vervangt,  indien  aanwezig,  altijd  de  gekozen  toestandsvergelijking  voor  het 
kruitgas  in  de  input-file,  die  dan  dus  een  dummy-rol  vervult. 

Dan  volgt  de  subroutine  TWPHIN.  Deze  bevat  de  beginwaarden  van  enkele  gas-  en 
kruitgrootheden,  namelijk  de  snelheid  (UXI  en  2,  UYl  en  2),  dichtheid  (RHOI  en  2), 
volume  fractie  (VOLFRl  en  2),  massa  (SAMI  en  2),  interne  energie  (SlEl  en  2,  om  de 
temperatuur  te  bepalen,  hier  is  chemische  energie  niet  in  opgenomen),  de  hoeveelheid 
reeds  verb»’and  kruit  (PARDEN)  en  een  parameter  die  aangeeft  of  het  kruit  al  brandt 
(IFBRN).  Deze  subroutine  is  verplicht  bij  Grain  burning,  en  de  INITIAL  kaarten  van  kruit 
en  kruitgas  zijn  hier  dan  ook  altijd  dummies. 

De  langste  subroutine  in  het  voorbeeld  is  EXSGVG,  waarin  de  oppervlakte-volumc 
verhouding  van  de  kruitkonrels  op  ieder  willekeurig  moment  wordt  berekend.  Standaard  is 
het  programma  namelijk  alleen  met  een  routine  uitgemst  die  dit  doet  voor  sferische 
korrels.  De  Otomelara  gebruikt  echter  cylindrisch  7-gats  kruit.  Deze  subroutine  moet 
worden  gedeclareerd  in  de  input-file  met  IFSGVG,  anders  wordt  de  standaard  routine 
(sferisch  kruit)  genomen.  In  de  routine  in  het  voorbeeld  (die  overigens  ook  voor  0-,  1-  of 
19-gats  kruit  te  gebruiken  is,  na  het  wijzigen  van  enkele  variabelen)  wordt  met  behulp  van 
de  subroutine  TABFRM  een  lijst  van  oppervlakte-volume  verhoudingen  berekend,  waarbij 
ook  rekening  wordt  gehouden  met  het  uiteen  vallen  van  de  korrel  in  slivers  na  het 
doorbranden  van  het  web  (subroutines  GENIS  en  GENOS).  Vervolgens  wordt  voor  een 
gegeven  kruitmassa  via  interpolatie  (subroutine  DLAGR)  het  bijbehorende  oppervlakte 
bepaald. 

Met  de  subroutine  EXEDIT  kan  additionele  informatie  worden  geprint.  In  het  voorbeeld 
wordt  de  maximale  druk  in  iedere  gridcel  en  in  het  hele  grid  gezocht  en  afgedrukt.  Deze 
subroutine  wordt  automatisch  na  iedere  cycle  opgeroepen. 

Ook  de  subroutine  EXBDFO  wordt  iedere  cycle  gebruikt.  Hierin  worden  extra  krachten  op 
het  projectiel,  ten  gevolge  van  het  loskomen  uit  de  huls,  het  insnijden  van  de  geleidings- 
band  in  de  trekken  en  velden  en  de  wrijving  in  de  loop,  bepaald.  Over  deze  krachten,  die 
overigens  vaak  als  tegendruk  worden  behandeld,  zijn  geen  bijzonder  nauwkeurige 
gegevens  bekend.  De  gebruikte  waarden  zijn  dan  ook  naar  analogie  met  andere  wapens 
gekozen. 

In  de  laatste  subroutine,  EXBUR,  wordt  een  alternatieve  burning  rate  gedefinieerd.  Zoals 
al  eerder  aan  bod  kwam,  kent  het  standaard  programma  alleen  de  wet  van  Vieille:  v=I3*p“, 
waarin  v  de  burning  rate  is  en  p  de  druk  (het  zou  beter  zijn  om  B*(p/po)“  te  gebruiken,  in 
verband  met  de  eenheid  van  B,  maar  helaas  doet  het  programma  dat  niet).  Wij  gebruiken 
hier  de  wet  van  Muraour-Aunis:  v=c+B*p.  Deze  subroutine  moet  in  de  input-file  gedecla¬ 
reerd  worden  met  IFBURN. 

Er  zijn  nog  veel  meet  subroutines  mogelijk  en  het  zal  duidelijk  zijn  dat  het  model  op  deze 
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manier  haast  onbeperkt  kan  worden  aangepast. 
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3  DE  GEGEVENS  VAN  HET  OTOMELARA  76/62  KANON 

De  Otomelara  16162,  van  Italiaanse  makelij,  is  een  veel  voorkome  ’  wapen  binneii  de 
NATO.  De  Koninklijke  Marine  gebruikt  ze  op  de  S-fregatten,  pri.  lair  tegen  lucht-  en 
zeedoelen.  In  het  kader  van  een  onderzoek  naar  mogelijke  mondingssnelheid  vergroting  is 
een  simulatie  model  opgesteld  met  behulp  van  het  in  de  vorige  paragrafen  besproken 
programma  PISCES  2DELK  Grain  burning.  De  input-file  en  de  gebruikte  externe 
subroutines  zijn  daar  ook  al  genoemd,  en  worden  gegeven  in  de  appeitdices  1  en  2. 


3.1  Gegevens  en  herkomst 


Wat  betreft  de  input  gegevens  is  het  niet  altijd  even  gemakkelijk  geweest  om  aan  de  juiste 
informatie  te  komen,  eenvoudigweg  omdat  dit  soort  problemen  meestal  veel  empirischer 
wordt  aangepakt  en  er  daarom  vaak  geen  fundamenteel  theoretische  grootheden  bekend 
zijn.  Bovendien  zijn  niet  alle  kruitlots  hetzelfde,  zeker  niet  als  ze  van  versehillende 
producenten  afkomstig  zijn. 

Alle  gegevens,  die  gebruikt  zijn,  omtrent  het  wapen  en  het  kruit  zijn  afkomstig  van  de 
Afdeling  Wapen-  en  Communicatie  Systemen/  Bureau  Geschut  en  Munitie  (DMKM/- 
WCS/GEMU)  op  de  Admiraliteit  in  Den  Haag  [3],  de  NATO  schootstafel  [4],  de  firma 
Muiden  Chemie  BV  en  de  firma  Diehl  GmbH  [5],  de  laatste  twee  via  bureau  GEMU.  Bij 
onderlinge  inconsequenties  is  contact  opgenomen  met  bureau  GEMU. 


Zoals  de  naam  al  zegt,  heeft  de  Otomelara  een  kaliber  van  76  mm,  en  is  de  loop  62 
kalibers  lang. 

In  de  huls  wordt  2,35  kg  kruit  gebruikt,  met  de  volgende  korrelgeometrie: 


lengte 

15.40  mm 

diameter 

6.78  mm 

diameter  gat 

0.04  mm 

binnen  web 

1.21  mm 

buiten  web 

1.22  mm 

aantal  gaten 

7 

ing  is: 

nitrocellulose 

83.4  % 

diphenylamine 

1.2  % 

dibutylphtalaat 

3.0  % 

dinitrotolueen 

9.0  % 

kaliumsulfaat 

2.4  % 

water 

0.7  % 

ethanol  (oplosmiddel) 

0.3  % 

ethylcentralite 

<0.05  % 

grafiet 

sporen 

Het  heeft  een  specifieke  massa 

gemeten  verbrandingswarmte 
ontstekings  temperatuur 


(stikstoggehalte  13.15  %) 


(niet  meegenomen  in  de  berekening) 
(niet  meegenomen  in  de  berekening) 
1570  kg/m’ 

3208  kJ/kg 
445  K  (172°C) 
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Deze  gegevens  zijn  afkomstig  van  Muiden  Chemic.  Diehl  geeft  iets  kleinere  afinctingen. 
Dit  heeft  enige  invloed  op  de  reactiesnelheid.  omdat  het  brandbaar  oppervlakte  andcrs  is. 
Als  ontstekingslading  wordt  157  g  kruit  type  TR54  en  18  g  zwart  buskruit  gebruikt.  Van 
het  TR54  kruit  hebben  wij  geen  gegevens,  daar  waar  dit  wordt  gebruikt  (in  de  subroutine 
EXFLW2,  waar  de  aandrijflading  als  een  inflow  van  materiaal  is  genuxielleerd)  is 
uitgegaan  van  zwart  buskruit. 

De  massa  van  de  granaat  is  6.3  kg. 

Wat  betreft  de  te  verwachten  mondings.snelheid  en  de  maximale  druk  wordt  op  dit  moment 
nog  exacte(re)  informatie  gezocht  door  bureau  GEMU,  de  waarden  liggen  hier  tussen  900 
m/s  en  925  m/s,  respectievelijk  300  MPa  en  360  MPa. 


3.2  Uitwerking  van  de  kruitgegevens 

Het  programma  PISCES  2DELK  vraagt  als  input  de  burning  rate.  Deze  moet  worden 
bepaald  met  de  bovenstaande  gegevens  in  combinatie  van  experimenteel  bepaalde 
verbrandingscurves. 

Eerst  moet  uit  de  samenstelling  van  het  kruit  een  aantal  andere  grootheden  worden 
berekend.  Dit  is  gedaan  met  het  thermochemische  programma  BLAKE  [7]  dat  op  basis 
van  de  chemische  reacties  en  hun  evenwichten  de  reactieproducten  (hoofdzakelijk 
kruitgassen)  bepaalt.  Voor  een  ladingsdichtheid  van  0.2  g/cm^  zijn  de  resultaten; 
adiabatische  vlamtemperatuur  2554  K 
maximale  theoretische  druk 
force 

molaire  massa  kruitgassen 
covolume  kruitgassen 
gamma  (Cp/cy) 
soortelijke  warmte  Cp 
De  kruitgassen  bestaan  hoofdzakelijk  uit 
CO  20.97  mol/kg 
Hj  7.o2  molAg 
HjO  6.7 1  mol/kg 
Nj  4.44  mol/kg 
COj  3.00  mol/kg 
KOH  0.23  mol/kg 
HjS  0. 1 1  mol/kg 

De  overige  kruitgassen  en  vaste  en  vloeibare  reactieproducten  (in  totaal  zijn  26  gassen  en 
drie  vaste  stoffen  als  reactieproducten  meegenomen  in  de  berekening)  beslaun  tezamen 
minder  dan  '/i%. 

De  burning  rate  wordt  afgeleid  uit  experimenten;  zogenaamde  ’gesloten  bom’  (closed 
vessel)  proeven.  Het  idee  is  eenvoudig:  in  een  afgesloten  bus  wordt  kruit  aangestoken  en 
de  verbranding  wordt  de  druk  gemeten  als  functie  van  de  tijd. 
iu..jt  de  NATO-standaard  gesloten  bom  van  de  Koninklijke  Marine  nog  niet  operationeel 
is,  werden  wij  in  de  gelegenheid  gesteld  de  experimenten  op  de  Koninklijke  Militaire 


232.6  MPa 

915.2  kJ/kg 

23.2  g/mol 
1.065  cmVg 
1.258 

42.1  J/Mol*K 
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School  in  Brussel  (Belgie)  te  doen  in  het  laboratorium  van  professor  R.  Meysinaiis. 

De  experimenten  zijn  uitgevoerd  met  kruittype  M6+2.  lotnummer  K-22/S/1 -BPD-79  in 
twee  verschillende  bommen:  het  volume  van  de  grooLste  is  6X2.3  cm',  dat  van  de  kleinste 
156.3  cm\  De  volume-oppervlakte  verhouding  van  beide  bommen  is  niet  hetzelfde.  de 
kleine  bom  heeft  relatief  een  grot  oppervlakte.  Het  kruit  werd  aangestoken  met  een 
gasmengsel  van  33%  methaan  en  6  zuurstof  op  een  totale  druk  van  0.3  MPa  (3bar).  dat 
tot  ontbranding  gebracht  dour  ne.  joorbranden  van  een  dun  stroomdraadje  (dat  dus  bij 
ieder  ’schot’  vemieuwd  moet  worden). 

De  volgende  experimenten  zijn  uitgevoerd: 


ladingsuichtheid  (g/cm^) 

grote  bom 

kleine  bom 

0.10 

M6-2 

M6-11,  M6-14 

0.15 

M6-3 

M6-12,  M6-15 

0.20 

M6-4 

M6-13,  M6-16 

0.25 

M6-5 

0.30 

M6-6 

De  maximale  druk,  gemeten  bij  M6-6,  was  373.7  MPa. 

Door  de  gemeten  druk-tijd  (p-t)  curves  is  vervolgens  een  polynoom  fit-curve  gelegd 
(maximaal  van  de  21e  graad,  meestal  ongevee,  lOe  graads).  Deze  fitcurves  zijn  daarna 
omgeschreven  naar  vivaciteits  (levendigheids)  curves,  volgens 

1 

dt  p-po 


Hierin  is  po  de  maximale  druk. 

In  figuren  1  en  2  is  deze  vivaciteitsfit  L  van  de  experimenten  uitgezet  tegen  p/p„. 

Vervolgens  moet  uit  de  meetwaarden  de  burning  rate  worden  bepaald.  Het  principe 
hiervan  is  als  volgt.  Er  vanuit  gaande  dat  tijdens  de  verbranding  het  kruit  wordt  omgezet 
tot  kruitgas  met  de  (constante)  adiabatische  vlamtemperatuur,  kunnen  we  stellen  dat  de 
toename  van  de  druk  evenredig  is  met  de  toename  van  de  hoeveelheid  kruitgas  en  dus  met 
de  afname  van  de  hoeveelheid  kruit.  Omdat  het  kruit  een  (nagenoeg)  constante  dichtheid 
heeft,  is  dus  het  volume  verbrand  kruit  per  tijdseenheid  evenredig  met  de  druktoename  per 
tijdseenheid.  Als  we  nu  het  oppervlakte  van  de  kruitkorrel  kennen  (hetgeen  voor  iedere 
geometrie  is  uit  te  rekenen)  weten  we  ook  de  inbranddiepte  per  tijdseenheid  of  burning 
rate,  omdat  de  inbranddiepte  gelijk  is  aan  de  volumeverandering  per  oppervlakte. 

In  werkelijkheid  moet  ook  nog  rekening  gehouden  worden  met  het  covolume  van  de 
kruitgassen,  de  onsteking  (beide  bekend)  en  warmteverliezen  naar  de  bom.  Van  dit  laatste 
kan  een  redelijke  schatting  worden  gemaakt  door  de  theoretische  en  de  gemeten  maximale 
druk  met  elkaar  te  vergelijken.  Het  verschil  hiertussen  is  (in  eerste  orde)  een  maat  voor 
verloren  warmte. 

Merk  op  dat  deze  berekening  werkt  met  de  p-t  curve  en  niet  met  de  vivaciteitscurve.  De 
vivaciteit  is  hier  slechts  gebruikt  om  de  metingen  grafi.sch  weer  te  geven,  en  de  gebruikte 
definitie  (soms  komen  varianten  voor)  doet  dus  in  de  berekening  van  de  burning  rate  niet 
ter  zake. 

Wij  hadden  twee  programma’s  tot  onze  beschikking  die  de  burning  rate  berekenden  uit  de 
meetwaarden,  namelijk  TBI12  [7|  en  BRL-BR  (8|.  In  figuur  3  zijn  de  resultaten  te  zien 
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Figuur  1:  Levendigheidscurves  van  Otomelara  kruit  gemeten  met 
de  grote  bom. 

van  berekeningen  met  BRL-BR.  Hierin  zijn  M6-2,  M6-11  en  M6-14,  met  een  ladingsdicht- 
heid  van  0.1  g/cm^,  weggelaten  voor  de  overzichtelijkheid  van  de  figuur.  Zij  gedroegen 
zich  geheel  zoals  de  andere  en  lagen  ’op’  de  andere  curves.  In  appendix  5  zijn  de  figuren 
1  t/m  3  in  kleur  en  op  A3-formaat  weergegeven,  zodat  ook  de  verschillende  metingen 
goed  te  onderscheiden  zijn. 

Burning  rates  berekend  met  TBI12  zijn  vrijwel  identiek. 

We  merken  op  dat  aan  het  einde  van  iedere  curve,  nadat  hij  is  begonnen  te  dalen,  een 
vreemde  ’bult’  te  zien  is.  Dit  is  een  artefact  van  de  berekening.  De  verklaring  is  dat  de 
kruitkorrels  op  dit  punt  uit  elkaar  beginnen  te  vallen,  alleen  niet  precies  allemaal  tegelij- 
kertijd.  Het  berekende  korreloppervlakte  gaat  dan  sterk  afwijken  van  de  werkelijke  waarde, 
zodat  de  daarmee  berekende  burning  rate  niet  meer  realistisch  is.  Een  sterke  aanwijzing 
hiervoor  is  dat  de  ’bult’  steeds  alleen  op  het  laatst  optreedt,  en  steeds  slechtst  een  procent 
of  15  in  beslag  neemt.  Dit  is  precies  het  volume  van  de  kruitkorrel  op  het  moment  dat  hij 
uiteen  gaat  vallen  ten  opzichte  van  de  oorspronkelijke  korrel  (13.4%). 

Ook  de  vorm  van  de  bult  is  hieruit  te  verklaren.  Daartoe  bekijken  we  figuur  4,  die  het 
oppervlakte  van  de  kruitkorrel  geeft  als  functie  van  het  volume,  zoals  dat  wordt  berekend 
door  de  beide  programma’s  TBI12  en  BRL-BR  (en  ook  in  de  exteme  subroutine  EXSGVG 
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Figuur  2:  Levendigheidscurves  van  Otomelara  kruit  gemeten  met 
de  kleine  bom. 


in  het  programma  PISCES  2DELK  Grain  burning). 

We  zien  een  zecr  duidelijke  knik  optreden  op  het  moment  dat  de  (7-gats)  korrel  in  slivers 
uit  elkaar  valt.  De  programma’s  gaan  er  vanuit  dat  dit  bij  alle  kruitkorrels  op  precies 
hetzelfde  moment  gebeurt.  Doordat  er  een  zekere  spreiding  is  in  de  afmetingen  van  de 
korrels  is  dat  onjuist.  Dus  voordat  het  programma  dat  verwacht,  zullen  sommige  korrels  al 
uiteen  vallen,  waardoor  hun  oppervlakte  enorm  afneemt  (zie  figuur  4).  Het  berekende 
totale  oppervlakte  is  dan  dus  groter  dan  het  werkelijke.  Omdat  de  burning  rate  het 
verbrande  volume  per  oppervlakte  (per  tijdseenheid)  is,  berekent  het  programma  dus  een  te 
kleine  burning  rate.  Dit  is  te  zien  in  figuur  3,  v66r  de  bult  daalt  de  burning  rate.  Aan  de 
andere  kant,  op  het  moment  dat  het  programma  verwacht  dat  alle  kruitkorrels  (tegelijker- 
tijd)  uil  elkaar  gevallen  zijn,  zullen  er  nog  korrels  intact  zijn.  Die  hebben  dus  een  groter 
oppervlakte  dan  het  programma  berekent,  en,  analoog  met  het  bovenstaande,  berekent  het 
programma  een  te  grote  burning  rate. 

Precies  dezelfde  redenering  is  van  toepassing  op  het  opbranden  van  de  slivers  kruit. 
Doordat  dit  ock  niet  voor  alle  slivers  simultaan  zal  gebeuren,  hoewel  het  programma  dat 
wel  verwacht,  daalt  de  burning  rate  uiteindelijk  en  houdt  hij  niet  abrupt  op. 
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Figuur  3 :  Burning  rate  curves  van  Otomelara  kruit . 

De  beide  programma’s  TBI12  en  BRL-BR  bieden  ook  de  mogelijkheid  om  aan  de  op 
bovenstaande  manier  gevonden  burning  rate  curves  een  verbrandingswet  te  fitten.  Helaas 
werken  beide  alleen  met  de  wet  van  Vieille:  v=B*p“,  waar  geen  constante  term  in 
mogelijk  is.  Wij  concluderen  echter  uit  de  grafiek  (figuur  3)  dat  een  rechte  lijn  die  niet 
door  de  oorsprong  loopt  geschikter  is,  hetgeen  neerkomt  op  de  wet  van  Muraour-Aunis: 
v=c+B*p.  We  bepalen  dan  c  en  B  als: 

c  (m/s)  6  (m/s*Pa) 

grote  bom  proeven  0.01  8.0*10  '° 

kleine  bom  proeven  0.01  9.0*10 '° 

We  zullen  in  de  simulatie  met  PISCES  2DELK  beide  waarden  gebruiken,  hoewel  de  grote 
bom  proeven  nauwkeuriger  zijn,  omdat  ze  in  veel  grotere  drukgebieden  zijn  gemeten  en 
omdat  de  oppervlakte  volume  verhouding  kleiner  is  en  dus  ook  de  waimteverliezen  naar 
de  wand. 

In  de  Powder  Description  Sheet  van  Muiden  Chemie  is  ook  een  vivaciteitscurve  gegeven. 
Deze  heeft  waarden  die  iets  lager  liggen  dan  de  door  ons  zelf  gemeten  curves,  namelijk 
rond  0,lbar'‘s  '.  Het  gaat  hier  echter  om  een  slecht  leesbare  kopie.  Nadere  informatie  was 
niet  beschikbaar.  Gewenste  informatie  is:  de  ’debulleting’  kracht  (om  granaat  uit  huls  los 
te  trekken),  de  dimensies  en  het  gebruikte  materiaal  van  de  geleidingsband,  een  diagram 
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Figuur  4:  Oppervlakte  van  een  Otomelara  kruitkorrel  als  func- 
tie  van  het  volume. 

van  de  aandrijvende  kracht  en  de  wrijvingskracht  in  de  loop,  de  vivaciteits  curve  van  het 
kruit  en  een  diagram  van  de  gasdruk  (versus  tijd  en  afgelegde  weg  in  de  loop)  met  de 
maximum  druk. 


3.3  Grootheden  die  onvoldoende  nauwkeurig  bekend  zijn 

Behalve  grootheden  waarover  nog  niet  voldoende  nauwkeurigheid  bestaat,  zoals  de  lostrek- 
kracht,  de  wrijving  in  de  loop,  de  insnijkracht,  de  maximale  druk,  de  mondingssnelheid,  de 
overeenkomst  tussen  gemeten  en  opgegeven  kruitgegevens,  onstekingslading,  etcetera, 
maar  die  met  behulp  van  bureau  GEMU  bekend  moeten  worden,  zijn  er  ook  nog  gegevens 
die  binnen  de  Koninklijke  Marine  of  bij  de  leverancier  onbekend  zijn.  Dit  zijn  met  name 
de  waarden  die  in  de  input  van  PISCES  2DELK  onder  GRNMISC  (zie  paragraaf  2.3) 
vermeld  moeten  worden.  Het  is  gelukt  om  hier  enige  fundamentele  literatuur  over  te 
verzamelen  [10],  waar  echter  wegens  onderlinge  inconsistentie  (en  tijdgebrek  om  daar  diep 
op  in  te  gaan)  nog  niet  genoeg  gebruik  van  is  gemaakt.  De  zaak  is  daarom  pragmatischer 
aangepakt,  en  van  de  genoemde  grootheden  is  bekeken  in  welke  mate  zij  de  uitkomsten 
van  de  berekening  bei'nvloeden.  Gelukkig  bleek  geen  van  de  onbekende  gegevens  van 
dramatische  invioed  op  de  resultaten.  Dit  zal  in  het  volgende  hoofdstuk  worden  beschre- 
ven. 
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SIMULATIE  VAN  EEN  SCHOT  IN  DE  OTOMELARA  76/62 


In  dit  hoofdstuk  worden  de  resultaten  van  de  simulatie  met  PISCES  2DELK  Grain  burning 
van  de  Otomelara  76/62  besproken.  Het  gebruikte  model  en  de  gegevens  zijn  in  de  vorige 
bladzijden  behandeld  en  terug  te  vinden  in  de  input-file  en  de  externe  subroutines, 
respectievelijk  appendix  1  en  2. 

De  simulatie  is  uitgevoerd  met  twee  verschillende  grids,  een  grove  en  een  fijnere.  Het 
grove  bestaat  uit  een  enkele  rij  van  50  cellen  (zie  figuur  5),  het  fijne  uit  6  rijen  van  99 


cellen  (figuur  6). 


Figuur  5:  Het  grove  Otomelara  grid,lL50. 


Figuur  6:  Het  fijne  Otomelara  grid,  6L100. 


Het  fijne  grid  is  het  grid  uit  appendix  1.  In  beide  gevallen  is  overigens  een  extra  rij 
gebruikt  om  de  vlampijp  met  de  ontstekingslading  te  modelleren.  Deze  wordt  verder  niet 
gebruikt.  We  zullen  de  berekeningen  met  beide  grids  aanduiden  met  1L50  voor  het  grove 
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grid  en  6L100  voor  het  fijne. 

In  de  simulatie  wordt  gewerkt  met  de  burning  rate  die  is  berekend  met  de  gegevens  uit  de 
grote  bom  (paragraaf  3.2):  v=c+B*p  met  c=0.0l  m/s  en  B=8.0*l()  ‘”nVs*Pa. 


4.1  Resultaten 

De  resultaten  zijn  als  volgt: 


grid 

duur 

aantal 

mondings- 

max.  druk 

tijdstip(ms) 

hoeveelheid 

(ms) 

cycles 

snelheid(m/s) 

(MPa) 

max  .druk 

kruitgas(kg) 

1L50 

9.46 

1119 

950 

475 

4.84 

2.496 

6L100 

9.47 

7479 

957 

509 

4.86 

2.515 

Merk  op  dat  inderdaad  de  gemiddelde  tijdstap  (duur  gedeeld  door  aantal  cycles)  van  het 
grove  grid  veel  groter  is  dan  die  van  het  fijne,  zoals  al  in  paragraaf  2. 1  werd  opgemerkt. 

In  beide  gevallen  treedt  de  maximale  druk  op  in  het  begin  van  de  huls,  bij  het  grove  grid 
ongeveer  15  cm  van  de  wand,  bij  het  fijne  grid  helemaal  bij  het  begin  direkt  boven  de 
vlampijp. 

De  mondingssnelheid  van  het  fijne  grid  is  iets  (0.7%)  groter  dan  die  van  het  grove  grid. 
Dit  kleine  verschil  komt  bijna  precies  overeen  met  het  verschil  in  de  totale  hoeveelheid 
kruitgas  (die  een  superpositie  is  van  de  verbrande  voortdrijvende  lading  en  de  uit  de 
vlampijp  naar  binnen  gestroomde  ontstekingsgassen). 

De  maximale  druk  geeft  een  veel  groter  verschil  (ruim  7%)  tussen  beide  grids  te  zien;  de 
maximale  druk  is  een  stuk  groter  bij  het  fijne  grid.  Dit  is  begrijpelijk.  Bij  het  fijne  grid 
wordt  namelijk  veel  meer  van  het  golfkaraktcr  van  de  drukoplx)uw  behoudcn.  Bij  het 
grove  grid  worden  de  extremen  veel  meer  uitgemiddeld.  Dit  is  goed  te  zien  in  de  figuren 
17  t/m  22  die  enkele  drukprofielen  gesimuleerd  met  de  2  verschillende  grids  tonen. 

Het  zou  interessant  zijn  om  in  de  lengte  richting  met  een  nog  fijner  grid  te  werken.  Het 
programma  kent  echter  een  maximum  van  KX)  gridcellen. 

In  de  volgende  bladzijden  (figuren  7  t/m  25)  worden  de  plots  van  beide  grids  met  diverse 
variabelen  als  functie  van  tijd  en  plaats  in  de  loop  steeds  naast  elkaar  gezet. 
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V(M/5)  VAN  HET  PROJECT I EL 


OTOMELORR  1L50 


COLUMN  I  RON  I  3U6GRI0  2 


0.0 

.002 

.  .  .  . 

.  .  .  . 

'  '  '  '  .'OM . 

TIJDtS) 

TIME 

Figuur  7:  De  snelheid  van  het  projectiel  in  de  loop  als  een 
functie  van  de  tijd  (berekend  met  1L50) . 
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Figuur  8:  De  snelheid  van  het  projectiel  in  de  loop 
als  een  functie  van  de  tijd  (berekend  met  6L100) . 
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V(M/5)  VAN  HET  PROJECTIEL 


OTOMELflRfl  1L50 


BOO. 


600. 


COLUMN  1  RON  1  3U6GRI0  2 


X 

PLflflTS  XCMI  VON  HET  PROJECTIEL  IN  DE  LOOP 


Figuur  9;  De  snelheid  van  het  projectiel  als  func 
tie  van  de  plaats  in  de  loop. 
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OTOMELflRfl  6L100 


COUIHN  I  ROW  1  SUBGRID  2 


X 

PLflflTS  XIHl  VAN  MET  PROJECTIEL 


Figuur  10:  Zie  figuur  9 


PLflfiTS  X(M) 


OTOMELflRfl  1L50 


COLUMN  I  ROM  I  SUBGRID  2 


Figuur  11:  De  plaats  van  het  projectiel  in  de  loop  als  een 
functie  van  de  tijd. 
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PLwrrs  x(M) 


OTOMELRRfl  6L100 


COLUMN  1  RON  1  SUBGRIO  2 


Figuur  12:  Zie  figuur  11. 
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OTOMELHRfl  6L100 


TIME  HISTC?!  FOR  MfTTERIRL  1 


2.0 


O 


TIME 

TlJOtSl 


Figuur  14:  Zie  figuur  13. 


OTOMELRRR  1L50 


PLPRTS  X(M)  IN  OE  LOOP 


uur  15;  Volume  fractie  kruii 
de  loop  op  elKe  milliseconds 


OTOMELflRfl  1L50 


Figuur  17:  Druk  als  een  functie  van  tijd  op  verschillende 
plaatsen  in  de  loop;  (1)  is  het  drukverloop  in  het  begin  van 
de  huls,  (5)  op  ongeveer  een  meter  van  het  einde  van  de  loop. 


0.6  '  .'ofi'  .00«  .008 

TlJOtSl 


Figuur  18:  Zie  figuur  17;  vanwege  het  andere  grid  zijn  de 
plaatsen  in  de  loop  waarvan  het  drukverloop  is  gegeven  iets 
anders  dan  in  figuur  17. 
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OTOMELRRfi  6L100 


COLUMN  1  RON  5  SUBGRID  1 


tu 
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TIME 
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Figuur  20:  Zie  figuur  19. 
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OTOMELflRfl  6L100 


in 


o 

■ 


E 


STRRflL  R(H)  IN  HET  BEGIN  VRN  OE  HULS 

K1 


Figuur  23:  De  druk  langs  de  straal  van  de  loop  (het  nulpunt 
stemt  overeen  met  de  bovenkant  in  het  begin  van  de  huls)  om  de 
miliseconde;  drukprof ielen  na  het  tijdstip  van  maximale  druk 
zijn  gestreept  weergegeven. 


De  figuren  23  en  24  tonen  radiale  drukprofielen  op  verschillende  plaatsen  in  de  huls.  Het 
nulpunt  komt  overeen  met  de  bovenkant  van  de  huls,  het  punt  op  ongeveer  4.5  cm  met  de 
bovenkant  van  de  vlampijp.  Deze  betreffen  alleen  het  fijne  grid  daar  het  grove  grid  slechts 
een  kolom  bevat. 

We  zien  dat  de  profielen  bijna  vlak  zijn  (alleen  boven  de  vlampijp  zit  een  "verstoring") 
wat  impliceert  dat  1  kolom  volstaat. 
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PtPfi)  VOOR  T=0+N«DT*  DT=1  MS 


OTOMELARfl  6L100 


OTOMELRRfl  6L100 


COLUMN  I  RON  2  aUBGRID  I 


u 


TlJOtSl 


Figuur  26:  Zie  figuur  25 
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Zoals  al  eerder  is  besproken  staan  er  in  de  input-file  enkele  onbekende  variabelen, 
namelijk  in  de  GRNMISC  kaarten.  We  hebben  daar  waarden  overgenomen  uit  een  studie 
van  MacNeal  Sehwendler  aan  het  Bofors  40/70  kanon  [11],  Om  toch  een  indruk  le  krijgen 
van  de  mate  waarop  deze  variabelen  het  resultaat  beinvloeden,  hebben  we  het  programma 
ook  laten  draaien  met  andere  getallen.  Het  resultaat  daarvan  is  als  volgt: 


grootheid 

variatie  grootheid 

variatie  mondingssnelheid 

prop,  const  interfacial  drag 

10'’  kleiner 

-3% 

critical  volume  fraction  solid 

-f-55% 

-1% 

laminar  viscosity  gas 

10^  groter 

+2% 

thermal  conductivity  gas 

10^  groter 

+27c 

thermal  conductivity  solid 

10‘  groter 

+  1% 

bulk  modulus  solid 

10  groter 

+  I7c 

specific  heat  solid 

10  groter 

-2% 

Uit  deze  tabel  blijkt  dat  deze  variabelen  de  berekening  nauwelijks  beinvloeden,  en  het  dus 
niet  onoverkomelijk  is  dat  zij  min  of  meer  onbekend  zijn.  Aan  de  andere  kant  kunnen  zij 
tezamen  een  fout  van  enkele  procenten  veroorzaken. 

Een  parameter  die  veel  beter  bekend  is  maar  desondanks  een  grote  fout  in  de  mondings- 
snelheid  en  met  name  de  maximale  druk  kan  veroorzaken  is  het  covolume.  Bij  een 
onnauwkeurigheid  in  deze  parameter  van  10%  varieert  de  mondingssnelheid  ongeveer  2% 
en  de  maximale  druk  10%. 


4.2  Discussie 

De  berekende  mondingssnelheid  en  maximale  druk  zijn  groter  dan  de  experimenteel 
gemeten  waardes.  Quantitatief  is  dit  precies  wat  men  verwacht:  er  is  namelijk  nog  geen 
rekening  gehouden  met  warmteverliezen  naar  de  kamer  en  de  loop  van  het  kanon.  Deze 
waimteverliezen  kunnen  oplopen  tot  een  procent  of  twintig  [12]  van  de  totale  beschikbare 
energie.  Hierdoor  daalt  zowel  de  mondingssnelheid  als  de  maximale  druk.  Het  is  boven- 
dien  aannemelijk  dat  de  maximale  druk  -aan  de  wand,  want  daar  zit  de  drukopnemer- 
meer  daalt  dan  de  projectielsnelheid.  Met  de  nieuwe  versie  van  Grainburning  die  ook  de 
warmteoverdracht  naar  de  wand  kan  simuleren  zal  dit  nader  onderzocht  worden. 

Verder  is  het  natuurlijk  mogelijk  dat  er  in  de  invoerparameters  en  dan  met  name  de  zo 
belangrijke  kruitverbrandingssnelheid  parameters  nog  te  grote  fouten  zitten. 

We  vermelden  verder  dat  een  vereenvoudigde  simulatie  van  de  gesloten  bom  proeven 
waarmee  de  burning  rate  werd  bepaald  druk-tijd  (of  levendigheids)  grafieken  opleverden 
die  zeer  goed  met  de  meetresultaten  overeenkomen.  In  deze  simulaties  werd  uiteraard  de 
experimenteel  bepaalde  burning-rate  gebruikt.  Het  is  bemoedigend  dat  been  en  terug 
rekenen  met  twee  programma’s  die  totaal  anders  van  aard  zijn  een  consistent  resultaat 
geeft. 
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APPENDIX  1 


De  input  file  (*.I2D) 


Het  grove  grid  1 L50 
ALLGRID  2  AXIAL  START 

HEADING  OTO  MELARA  25-2-91  7-GATSROUTINE  ENKEL  GRID  1L50 
COMMENT  *************************************** ********* 

COMMENT  *  * 

COMMENT  *  GRAIN  BURNING  IN  A  GUN  BARREL  SYSTEM  * 
COMMENT  *  ♦ 

COMMENT  ************************************************ 

COMMENT  1234567  1234567  1234567  1234567  1234567  1234567  1234567  1234567 
WRAPUP  l.OE-11  10000  900000  .2 

TSTEP  5.0E-7  0.67 

COMMENT - IDEAL  GAS  E.O.S.  FOR  GASPHASE  (MPN=1) — 

MATERIALBASIC  GAS  0.86728  0. 

MATERIALEOSA  GAS  GAMMA  1.258 

COMMENT . IDEAL  GAS  E.O.S.  FOR  GRAINS  (MPN=2)— - 

MATERIALBASIC  GRAIN  1570.  0. 

MATERIALEOSA  GRAIN  GAMMA  1.258 

COMMENT . INITIAL  VALUES - - 

COMMENT 

COMMENT . GAS - 

INITIAL  1  DENSITY  0.86728  SIE  4.1934E5XVEL  0.  YVEL  0. 

COMMENT . SOLID - 

INITIAL  2  DENSITY  1570.  SIE  5.8800E5XVEL  0.  YVEL  0. 

COMMENT . 

COMMENT 

COMMENT . RIGID - - 

INITIAL  4  XVEL  0.  YVEL  0.  BODYMASS6.30 
COMMENT 

COMMENT  — . - . 

POLYGON  SETUP  PROJECT  5 
COMMENT . - 

POLYGON  SETUP  FLOW  5  0.038  0.038 


POLYGON  POINTS 

FLOW 

1 

-.001 

-.001 

POLYGON  POINTS 

FLOW 

2 

.23 

-.001 

POLYGON  POINTS 

FLOW 

3 

.23 

.000001 

POLYGON  POINTS 

FLOW 

4 

.16 

.00901 

POLYGON  POINTS 

FLOW 

5 

-.001 

.00901 

COMMENT 


ALLEDIT  PRINT  NEVER 
ALLEDIT  ARCHIVE  NEVER 
ALLEDIT  INFSUM  NEVER 


ALLEDIT  EBDSUM  NEVER 
ALLEDIT  SUBSUM  NEVER 

ALLEDIT  LAGSUM  BYCYCLES  6000  WRAPUP  6000 
ALLEDIT  EULSUM  BYCYCLES  6000  WRAPUP  6000 
ALLEDIT  MATSUM  BYCYCLES  00  WRAPUP  25 
ALLEDIT  RESTART  BYCYCLES50000  WRAPUP  50000 

COMMENT . . — - - 

CUTOFF  5.E-4  4000.  l.E-14 

COMMENT - - - - 

COMMENT  . . INPUT  DATA  TWO  PHASE  FLOW- . 

COMMENT  . . . 

COMMENT  — -  (I  )  THERMAL  CONDUCTIVITY  GAS 
COMMENT  (2  )  SPECIFIC  HEAT  GAS 
COMMENT  (3  )  LAMINAR  VISCOSITY  GAS 

COMMENT  —  (4  )  PROPORTIONALITY  CONSTANT  INTERFACIAL  DRAG 

COMMENT  --  (5  )  THERMAL  CONDUCTIVITY  SOLID 

COMMENT  (6  )  SPECIFIC  HEAT  SOLID 

COMMENT  —  (7  )  BULKMODULUS  SOLID 

COMMENT  —  (8  )  CRITICAL  VOLUME  FRACTION  SOLID 

COMMENT  —  (9  )  INITIAL  PARTICLE  DIAMETER 

COMMENT  —  (10)  COMBUSTION  ENERGY  SOLID 

COMMENT  —  (11)  PARTICLE  BURNING  RATE  PROPORTIONALITY  CONSTANT 

COMMENT  — -  (12)  PARTICLE  BURNING  RATE  INDEX 

COMMENT  — -  (13)  IGNITION  TEMPERATURE 

COMMENT  —  (14)  QUENCHING  TEMPERATURE 

COMMENT . 

COMMENT  1234567  1234567  1234567  1234567  1234567  1234567  1234567  1234567 
GRNMISC  1  0.2291  2  1823.6  3  8.98 lE-5  4  1.75 

GRNMISC  5  0.2  6  2000.  7  1.3781E9  8  0.55 

GRNMISC  9  0.002  10  3.21E6  11  5.8E-8  12  0.775 

GRNMISC  13  445.  14  294.  15  10. 

COMMENT . - - - - - 

EXTRAS  1  0.5791  2  1.258  3  34.E6  4  2.E6 

COMMENT  EXTRA(l)  IS  THE  INITIAL  VOLFRl 
COMMENT  EXTRA(2)  IS  THE  GAMMA 

COMMENT  EXTRA(3)  AND  (4)  ARE  THE  PF2  AND  PF4  FRICTION  PRESSURES 

COMMENT  OF  THE  PROJECTILE  IN  EXBDFO 

ENDECK 

SUBGRID  EULER  TUBE  3  51 

COMMENT . . . 

COMMENT . STANDARD  MATERIAL  MODELS  TWO  PHASE  FLOWS 

COMMENT . . . 

IFTWPH 

IFBURN 

IFSGVG 

COMMENT  IFTEMP  *TSTNER* 

COMMENT . - . - 

ZONING  POINTS  1  1  0.01350  0.05000 


ZONING  POINTS  2  1  0.01350  0.00900 

ZONING  POINTS  1  3  0.16000  0.04725 

ZONING  POINTS  2  3  0.16000  0.00900 

ZONING  POINTS  1  4  0.23000  0.04600 

ZONING  POINTS  2  4  0.23000  0.000001 

ZONING  SEGMENT  ILINE  2  1  3 

ZONING  SEGMENT  ILINE  1  I  3 

COMMENT - - - 

ZONING  POINTS  1  9  0.62840  0.03800 

ZONING  POINTS  2  9  0.62840  0.000001 

ZONING  SEGMENT  ILINE  1  4  9 

ZONING  SEGMENT  ILINE  2  4  9 

COMMENT . . 

ZONING  POINTS  1  50  4.7000  0.03800 

ZONING  POINTS  2  50  4.7000  0.000001 

ZONING  SEGMENT  ILINE  1  9  50  1.033 

ZONING  SEGMENT  ILINE  2  9  50  1.033 

COMMENT - - - 

ZONING  POINTS  3  1  0.0135  0.00000 

ZONING  POINTS  3  4  0.23000  0.00000 

ZONING  POINTS  3  9  0.62840  0.00000 

ZONING  POINTS  3  50  4.7000  0.00000 

ZONING  SEGMENT  ILINE  3  1  4 

ZONING  SEGMENT  ILINE  3  4  9 

ZONING  SEGMENT  ILINE  3  9  50  1.033 

COMMENT - 

ZONING  POINTS  1  51  5.7000  0.03800 

ZONING  POINTS  2  51  5.7000  0.000001 

ZONING  POINTS  3  51  5.7000  0.00000 

COMMENT . 


REGION 

2 

GRAIN  RECTANG  -1. 

8. 

0. 

0.20 

1 

REGION 

1 

GAS  RECTANG  -1. 

8. 

-0.01 

0.01 

2 

REGION 

1 

GAS  RECTANG  -1. 

8. 

0.03 

0.05 

2 

REGION  1 
COMMENT  -- 

GAS  RECTANG  -1. 

8. 

0.07 

0.09 

2 

REGION 

UNUSED  INDEX  2 

3 

1 

51 

5 

COMMENT . 

BOUNDARYFLOW  ILINE  2  1  4  TWOPHA  EXFLOW 

BOUNDARYWRAPUP  JLINE  51  1  2  XVEL  0.  0. 

COMMENT . - . . 

EDIT  PRINT  BYTIMES  0  WRAPUP  0.5E-3  EULPRT 

EDLIST  VNAMES  PRESSUREVOLFRl  VOLFR2  PARDEN  EXEULl 
EULPRT 

EDLIST  VNAMES  SIE  DENSITY  EULPRT 

COMMENT .  - 

COMMENT  GLOBAL  ARCHIVE 

COMMENT . 

EDIT  ARCHIVE  BYTIMES  0  WRAPUP  0.1  E-3  EULSAV2 
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EDLIST  VNAMES 
EDLIST  VNAMES 

COMMENT . 

COMMENT . 

COMMENT 

COMMENT . 

EDIT  ARCHIVE  BYTIMES  0 


PRESSUREVOLFR2  X 
TEMPI  TEMP2 


GLOBAL  ARCHIVE 


EULSAV2 

EULSAV2 


WRAPUP  0.1  E-3 


EULSAV3 


EDLIST  COLS  1 

EULSAV3 

EDLIST  ROWS  1  2 

3  4 

5 

EULSAV3 

EDLIST  VNAMES  EXEULl 

COMMENT . . . . 

ENDECK 

SUBGRID  RIGID  GESCHUT  1 

5 

EULSAV3 

ZONING  POINTS  1  1 

0.5880 

0.0590 

ZONING  POINTS  1  2 

0.5880 

-0.001 

ZONING  POINTS  1  3 

6.0000 

-0.001 

ZONING  POINTS  1  4 

6.0000 

0.0590 

ZONING  POINTS  1  5 

COMMENT . 

0.5880 

0.0590 

BOUNDARYINTERPOLILINE 

1  1 

5  PROJECT  1 

COMMENT . 

REGION  4  INDEX 

COMMENT . 

1  1 

1  5 

EDIT  PRINT  BYTIMES  0 

WRAPUP  0.5E-3 

RIGS  A  VI 

EDLIST  COLS  1 

RIGS  AVI 

EDLIST  ROWS  1 

RIGS  AVI 

EDLIST  VNAMES  X  XVEL 

RIGSAVl 

COMMENT  — . 

EDIT  ARCHIVE  BYTIMES  0 

WRAPUP  0.1  E-3 

RIGSAV2 

EDLIST  COLS  1 

RIGSAV2 

EDLIST  ROWS  1 

RIGSAV2 

EDLIST  VNAMES 

ENDECK 

ENDINPUT 


XVEL 


RIGSAV2 


Het  fijne  grid  6L100 


ALLGRID  2  AXIAL  START 

HEADING  OTO  MELARA  29-1-91  7-GATSROUTINE  FUN  GRID  6L100 
COMMENT  ************************************************ 

COMMENT  *  * 

COMMENT  *  GRAIN  BURNING  IN  A  GUN  BARREL  SYSTEM  * 
COMMENT  *  * 

COMMENT  ************************************************ 

COMMENT  1234567  1234567  1234567  1234567  1234567  1234567  1234567  1234567 
WRAPUP  LOE-9  900000  900000  .2 

TSTEP  5.0E-7  0.67 

COMMENT . - . IDEAL  GAS  E.O.S.  FOR  GASPHASE  (MPN=1)-— 
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MATERIALBASIC  GAS  0.86728  0. 

MATERIALEOSA  GAS  GAMMA  1.258 

COMMENT  . . -—IDEAL  GAS  E.O.S.  FOR  GRAINS  (MPN=2) . 

MATERIALBASIC  GRAIN  1570.  0. 

MATERIALEOSA  GRAIN  GAMMA  1.258 

COMMENT . . INITIAL  VALUES - 

COMMENT 

COMMENT . GAS - - 

INITIAL  1  DENSITY  0.86728  SIE  4.1934E5XVEL  0.  YVEL  0. 

COMMENT  — . —SOLID - 

INITIAL  2  DENSITY  1570.  SIE  5.8800E5XVEL  0.  YVEL  0. 

COMMENT - - - 

COMMENT 

COMMENT . — RIGID - 

INITIAL  4  XVEL  0.  YVEL  0.  BODYMASS6.30 
COMMENT 

COMMENT . . 

POLYGON  SETUP  PROJECT  5 

COMMENT . 

POLYGON  SETUP  FLOW  5 
POLYGON  POINTS  FLOW  1 
POLYGON  POINTS  FLOW  2 
POLYGON  POINTS  FLOW  3 
POLYGON  POINTS  FLOW  4 
POLYGON  POINTS  FLOW  5 

COMMENT - - 

ALLEDIT  PRINT  NEVER 
ALLEDIT  ARCHIVE  NEVER 
ALLEDIT  INFSUM  NEVER 
ALLEDIT  EBDSUM  NEVER 
ALLEDIT  SUBSUM  NEVER 

ALLEDIT  LAGSUM  BYCYCLES  6000  WRAPUP  6000 
ALLEDIT  RESTART  BYCYCLES50000  WRAPUP  50000 


COMMENT  - — . 

CUTOFF  5.E-4  4000.  l.E-14 

COMMENT . 

COMMENT  - . INPUT  DATA  TWO  PHASE  FLOW — . 

COMMENT - 

COMMENT  —  (1  )  THERMAL  CONDUCTIVITY  GAS 
COMMENT  —  (2  )  SPECIFIC  HEAT  GAS 
COMMENT  —  (3  )  LAMINAR  VISCOSITY  GAS 

COMMENT  —  (4  )  PROPORTIONALITY  CONSTANT  INTERFACIAL  DRAG 


COMMENT  —  (5  )  THERMAL  CONDUCTIVITY  SOLID 
COMMENT  —  (6  )  SPECIFIC  HEAT  SOLID 
COMMENT  — -  (7  )  BULKMODULUS  SOLID 

COMMENT  —  (8  )  CRITICAL  VOLUME  FRACTION  SOLID 
COMMENT  —  (9  )  INITIAL  PARTICLE  DIAMETER 
COMMENT  —  (10)  COMBUSTION  ENERGY  SOLID 


0.038  0.038 

-.001  -.001 
.23  -.001 

.23  .000001 

.22450  .00901 
-.001  .00901 
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COMMENT  —  (11)  PARTICLE  BURNING  RATE  PROPORTIONALITY  CONSTANT 

COMMENT  —  (12)  PARTICLE  BURNING  RATE  INDEX 

COMMENT  (13)  IGNITION  TEMPERATURE 

COMMENT  —  (14)  QUENCHING  TEMPERATURE 

COMMENT - - - 

COMMENT  1234567  1234567  1234567  1234567  1234567  1234567  1234567  1234567 


GRNMISC  1 

0.2291  2 

1823.6  3 

8.98 lE-5  4 

1.75 

GRNMISC  5 

0.2  6 

2000.  7 

1.3781E9  8 

0.55 

GRNMISC  9 

0.002  10 

3.21E6  11 

5.8E-8  12 

0.775 

GRNMISC  13 

445.  14 

294.  15 

10. 

COMMENT  — - 
EXTRAS  1 

0.5809  2 

1.258  3 

34.E6  4 

2.E6 

COMMENT  EXTRAS(l)  IS  THE  INITIAL  VOLFR(l) 

COMMENT  EXTRAS(2)  IS  THE  GAMMA 

COMMENT  EXTRAS(3)  AND  (4)  ARE  THE  PF2  AND  PF4  PRESSURES 

COMMENT  OF  THE  PROJECTILE  IN  EXBDFO 

ENDECK 

SUBGRID  EULER  TUBE  8  100 

COMMENT  - . — 

COMMENT . STANDARD  MATERIAL  MODELS  TWO  PHASE  FLOWS 

COMMENT - 

IFTWPH 

IFSGVG 

IFBURN 

COMMENT  IFDRAG 
COMMENT  IFTEMP  ♦TSTNER* 

COMMENT . 

ZONING  POINTS  1  1  0.01180  0.04390 

ZONING  POINTS  2  1  0.00310  0.03760 

ZONING  POINTS  3  1  0.00000  0.02825 

ZONING  POINTS  4  1  0.00350  0.01880 

ZONING  POINTS  5  1  0.01060  0.01530 

ZONING  POINTS  6  1  0.01350  0.0 1 190 

ZONING  POINTS  7  1  0.01350  0.00900 

ZONING  POINTS  1  3  0.03410  0.04780 

ZONING  POINTS  1  5  0.05690  0.04950 

ZONING  POINTS  1  15  0.22450  0.04600 

ZONING  POINTS  7  15  0.22450  0.00900 

ZONING  SEGMENT  ILINE  7  1  15 

ZONING  SEGMENT  ILINE  1  1  3 

ZONING  SEGMENT  ILINE  1  3  5 

ZONING  SEGMENT  ILINE  1  5  15 

ZONING  SEGMENT  JLINE  15  1  7 

ZONING  BLOCK  AUTO  1  7  1  15 

COMMENT . - . - 

ZONING  POINTS  7  16  0.23000  0.000001 

ZONING  POINTS  1  16  0.23000  0.04600 

ZONING  SEGMENT  JLINE  16  1  7 


COMMENT  - . - . . 

ZONING  POINTS  1  33  0.62840  0.03800 

ZONING  POINTS  7  33  0.62840  0.000001 

ZONING  SEGMENT  ILINE  1  16  33 

ZONING  SEGMENT  ILINE  7  16  33 

ZONING  SEGMENT  JLINE  33  1  7 

ZONING  BLOCK  IDRAW  1  7  16  33 

COMMENT - 

ZONING  POINTS  1  99  4.7000  0.03800 

ZONING  POINTS  7  99  4.7000  0.000001 

ZONING  SEGMENT  JLINE  99  1  7 

ZONING  SEGMENT  ILINE  1  33  99  1.033 

ZONING  SEGMENT  ILINE  7  33  99  1.033 

ZONING  BLOCK  IDRAW  1  7  33  99 

COMMENT - 

ZONING  POINTS  8  1  0.0135  0.00000 

ZONING  POINTS  8  15  0.22450  0.00000 

ZONING  POINTS  8  16  0.23000  0.00000 

ZONING  POINTS  8  33  0.62840  0.00000 

ZONING  POINTS  8  99  4.7000  0.00000 

ZONING  SEGMENT  ILINE  8  1  15 

ZONING  SEGMENT  ILINE  8  16  33 

ZONING  SEGMENT  ILINE  8  33  99  1.033 

COMMENT - 

ZONING  POINTS  1  100  ^^7000  0.03800 

ZONING  POINTS  7  100  5.7000  0.000001 

ZONING  POINTS  8  100  5.7000  0.00000 

ZONING  SEGMENT  JLINE  100  1  7 


COMMENT 


REGION 

2 

GRAIN 

RECTANG  -1. 

8. 

0. 

0.20 

1 

REGION 

1 

GAS 

RECTANG-I. 

8. 

-0.01 

0.01 

2 

REGION 

1 

GAS 

RECTANG -1. 

8. 

0.03 

0.05 

2 

REGION  1 
COMMENT  - 

GAS 

RECTANG  -1. 

8. 

0.07 

0.09 

2 

REGION 

UNUSED  INDEX  7 

8 

1 

100 

5 

COMMENT - 

BOUNDARYFLOW  ILINE  7  1  15  TWOPHA  EXFLOW 

BOUNDARYWRAPUP  JLINE  100  1  6  XVEL  0.  0. 

COMMENT - 

EDIT  PRINT  BYTIMES  0  WRAPUP  0.5E-3  EULPRT 

EDLIST  COLS  3  6  EULPRT 

EDLIST  VNAMES  PRESSUREVOLFRl  VOLFR2  PARDEN  EXEULl 
EULPRT 

EDLIST  VNAMES  SIE  DENSITY  EULPRT 

COMMENT - 

COMMENT  GLOBAL  ARCHIVE 

COMMENT - 

EDIT  ARCHIVE  BYTIMES  0  WRAPUP  0.1  E-3  EULSAV2 
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GLOBAL  ARCHIVE 


EDLIST  VNAMES  PRESSUREVOLFR2  X 

COMMENT - 

COMMENT 

COMMENT - 

EDIT  ARCHIVE  BYTIMES  0 
EDUST  COLS  5 
EDLIST  VNAMES  TEMPI  TEMP2 

COMMENT - 

COMMENT - 

COMMENT  GLOBAL  ARCHIVE 

COMMENT - 

EDIT  ARCHIVE  BYTIMES  0 


WRAPUP  O.lE-3 


EULSAV2 


EULSAV4 

EULSAV4 

EULSAV4 


WRAPUP  O.lE-3 


EULSAV3 


EDLIST  COLS  1 

EULSAV3 

EDLIST  ROWS  1  2 

3  4 

5 

EULSAV3 

EDLIST  VNAMES  EXEULl 

EULSAV3 

COMMENT - - 

— 

ENDECK 

SUBGRID  RIGID  GESCHUT  1 

5 

ZONING  POINTS  1  1 

0.5880 

0.0590 

ZONING  POINTS  1  2 

0.5880 

-0.001 

ZONING  POINTS  1  3 

6.0000 

-0.001 

ZONING  POINTS  1  4 

6.0000 

0.0590 

ZONING  POINTS  1  5 

0.5880 

0.0590 

COMMENT - 

BOUNDARYINTERPOLILINE 

1  1 

5  PROJECT  1 

COMMENT  - . 

REGION  4  INDEX 

1  1 

1  5 

COMMENT . 

EDIT  PRINT  BYTIMES  0 

WRAPUP  0.5E-3 

RIGS  AVI 

EDLIST  COLS  1 

RIGSAVl 

EDLIST  ROWS  1 

RIGS  AVI 

EDLIST  VNAMES  X  XVEL 

RIGSAVl 

COMMENT . . 

EDIT  ARCHIVE  BYTIMES  0 

WRAPUP  0.1  E-3 

RIGSAV2 

EDLIST  COLS  1 

RIGSAV2 

EDLIST  ROWS  1 

RIGSAV2 

EDLIST  VNAMES 

ENDECK 

ENDINPUT 


XVEL 


RIGSAV2 
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APPENDIX  2 


De  exteme  subroutines 


PROGRAM  ELK2D  ELK2D  5 

C  *(INPUT,OUTPUT,TAPE5=INPUT,TAPE6  =  OUTPUT,TA- 

PE1,TAPE2,TAPE3,TAPE4,  ELK2D  11 

C  *  TAPE7.TAPE10,TAPE87,TAPE98;iAPE99)  ELK2D 

12 

C  ELK2D  14 


ELK2D  15 

C  ELK2D  16 

C  THE  MAIN  PROGRAM  IS  WRITTEN  BY  THE  USER  FOR  THE  FOLLOWING 
REASONS.  ELK2D  17 

C  THE  FILES  NEEDED  TO  EXECUTE  HIS  PROBLEM  ARE  ON  THE  PROGRAM 
CARD.  ELK2D  18 

C  THE  AMOUNT  OF  SMALL  CORE  MEMORY  AVAILABLE  IS  LENSCM. 

ELK2D  19 

C  THE  AMOUNT  OF  LARGE  CORE  MEMORY  AVAILABLE  IS  LENLCM. 

ELK2D  20 

C  THIS  IS  THE  DEFAULT  MAIN  PROGRAM  IF  THE  USER  DOES  NOT  SUPPLY 
ONE.  ELK2D  21 

C  ELK2D  22 

^  :|e  :tc  llclk  ifl  ifi  *  *  *  4<  *  ifc  4c  :ii  :|t  *  :k  l|c  * :)!  ikiiciK  4!  i|t  4c  ]|c  4c « 4c  4: 3<c  !(i  *  :|c  ;ic  *  *  *  :tc  III  :|c  *  )ii  :tc  4i  *  *  ♦ 


ELK2D  23 
C 
C 

COMMON  /CORE/  LIMCOR,  CORE(7000) 

26 

COMMON  /MEMORY/  LENSCM,  LENLCM.  GRID 
ELK2D  27 

C*IF  BEGIN  DEF,.NOT.(ECS.OR.LCM) 

28 

DIMENSION  GRID(30  000) 

31 

C*IF  END  DEF,.NOT.(ECS.OR.LCM) 

33 

LENSCM  =  0 
LENLCM  =  10  000  000 
C*IF  BEGIN  DEF,.NOT.(ECS.OR.LCM) 

36 


ELK2D  24 

ELK2D  25 

ELK2D 


ELK2D 

ELK2D 

ELK2D 

ELK2D  34 
ELK2D  35 
ELK2D 


LENSCM  =  30  000 
LENLCM  =  0 

C*IF  END  DEF,.NOT.(ECS.OR.LCM) 
42 
C 


ELK2D  39 
ELK2D  40 
ELK2D 

ELK2D  43 
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LIMCOR  =  7000 
CALL  CODE 
C 

C  END  OF  PROGRAM  MAIN 
47 

STOP 

END 

SUBROUTINE  EXFLW2  (UXEX,UYEX,ROEX,EINEX, 


ELK2D  44 

ELK2D  45 

ELK2D  46 

ELK2D 

ELK2D  48 

ELK2D  49 


1  UX1EX,UX2EX,UY1EX,UY2EX.R01EX,R02EX. 


C 

C 


2  SIE1EX,SIE2EX,VFR1EX.VFR2EX.PDENEX, 

3  PEX,IFUOPX^BAR,YBAR.TiM,I,J.M,MASMOM, 

4  NAMARY, 

5  VOLTOT,FCVNEW,AREA,FCSNEW. 

6  SXUNCSYUNC ) 

DIMENSION  NAMARY(2) 


><:****************************** 


c 

C  EXFLW2  IS  THE  USER  SUPPLIED  EXTRA  EULER  FLOW  BOUNDARY 
CONDITION 

C  FOR  TWO  PHASE  FLOW 

C 

C  QUANTITIES  TO  BE  COMPUTED  ARE 
C  PEX  -  PRESSURE 

C  UXIEX  -  X  COMPONENT  OF  FLUID(GAS)  VELOCITY 

C  UX2EX  -  X  COMPONENT  OF  SOLID  VELOCITY 

C  UYIEX  -  Y  COMPONENT  OF  FLUID(GAS)  VELOCITY 


C  UY2EX  -  Y  COMPONENT  OF  SOLID  VELOCITY 


C  ROIEX  -  FLUID(GAS)  MASS  DENSITY 

C  R02EX  -SOLID  MASS  DENSITY 

C  SIEIEX  -  FLUID(GAS)  SPECIRC  INTERNAL  ENERGY 


C  SIE2EX  -  SOLID  SPECIFIC  INTERNAL  ENERGY 

C  VFRIEX  -  FLUID(GAS)  VOLUME  FRACTION 
C  VFR2EX  -  SOLID  VOLUME  FRACTION 

C 

C  SWITCHES  (ALWAYS)  - 

C  IFUOPX  -  SWITCH,  INDICATING  TYPE  OF  BOUNDARY 
C  =0  -  ALL  FLOW  VARIABLES  ARE  DEFINED 
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c 


=1  -  PEX  .  VFRIEB  AND  VFR2EB  ARE  NOT  DEFINED 


C  =2  -  UXEX  AND  UYEX  (OR  UX1EX.UX2EX.UY1EX  AND  UY2EX)  ARE 

NOT  DEFINED 
C 

C  MASMOM  -  SWITCH  .  INDICATING  COLUMN  BUFFER  NUMBER 
C  =1  -  COLUMN  BUFFER  NUMBER  IB 

C  =2  -  COLUMN  BUFFER  NUMBER  ID 

C 

C  INTUT  ARGUEMENTS  ARE  - 
C  XBAR  -  AVERAGE  X  OF  FLOW  SEGMENT 

C  YBAR  -  AVERAGE  Y  OF  FLOW  SEGMENT 

C  TIM  -  TIME 

C  I,J,M  -  COLUMN,  ROW,  AND  SUBGRED  NUMBER 
C 

Q  :|c  iciic  :|c  *  :|c  **  !|e  **  lie  ;ii Ik  l|c  :|c  Ifclfc*  l|c  ****  !|i  **  1)1  «***:(<  4:  **  41  4c  **  ****«*******%;!);  4c  *  ](c  *  li!  !|c  **  * 

c 

C  THIS  COMMON  BLOCK  ASSUMES  LIMFLO  =  6  AND  LIMWAL  =  5 
EULBND  2 

C  THE  DIMENSION  OF  THESE  ARRAYS  IS  LIMEB  =  LIMFLO  +  LIMWAL  +  1 
EULBND  3 

COMMON  /EULBND/  NTREBD,  RADOPT,  GAMOUT,  SSOUT,  NUMFLO, 
IFUOP(12),  EULBND  4 

1  UXEB(12),  PREB(12),  XIMEB(12),  TXMEB(12),  TMAEB(12),  EINEB(12), 

EULBND  5 

2  UYEB(12),  ROEB(12),  YIMEE(12),  TYMEB(12).  TENEB(I2),  WRKEB(12), 

EULBND  6 

3  UX1EB(12),UX2EB(12),  UY1EB(12),  UY2EB(I2),S1E1EB(12),SIE2EB(12), 
EULBNDEX  1 

4  RO 1  EB(12),R02EB(  1 2), VFR 1  EB(  1 2),VFR2EB(  12),PDENEB(  1 2), 

EULBNDEX  2 

5  NAMEB(2,12),  LIMFLO,  LIMWAL,  LIMEB 
EULBNDEX  3 

C  EULBND  8 

C 

COMMON  /NCYVAR/  NCYCLE,  DLTMIN,  NSEXIT,  NUMPLT,  TIMB,  DLTB, 
DLTH  $/NCYVAR  2 

A  ,  DDSTEP,  ITSTEP,  KTSTEP,  NSWRAP,  NUMWAG,  TIME,  DLTE,  DLTFE 
$/NCYVAR  3 

B  ,  SSSTEP,  JTSTEP,  MTSTEP  $/NCYVAR 

4 

C  $/NCYVAR  5 

C 

COMMON  /COLVAR/  NCOLC,  EINTC(5),  XMOMC(5),  lA,  ID, 

COLVAR  2 

1  NCOLA,  NCOLD,  EKINC(5),  YMOMC(5),  IB,  IE, 

COLVAR  3 
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2  NCOLB,  NCOLE,  EDISC(5),  ZMASC(5),  EHTRH(5),  IC 
COLVAR  4 

C  COLVAR  5 

C 

C  THIS  DIMENSION  STATEMENT  ASSUMES  LIMROW(2)=100  AND  NBCOLS(2- 
)=5  EULARY  2 

COMMON  /COLMEM/  LIMCOB,  EULARY 

3 

AX  (100,5),  Y  (100,5),  THICK  (100,5),  UX  (100,5).  EULARY 

4 

B  UY  (100,5),  XMOM  (100,5),  YMOM  (100,5),  VOLMAS(  100,5), 

EULARY  5 

C  VOID  (100,5),  RHOl  (100,5),  RH02  (100,5),  SAMI  (100,5), 

EULARY  6 

D  SAM2  (100,5),  SIEl  (100,5),  SIE2  (100,5),  P  (100,5),  EULARY 

7 

E  Q  (100,5),  SSPD  (100,5),  DIV  (100,5),  ERGl  (100,5),  EULA¬ 

RY  8 

F  ERG2  (100,5),  NBTYPE(  100,5),  NGATES(  100,5),  FCOVV  (100,5), 

EULARY  9 

G  FCOVSI(  100,5),  FCOVSJ(  100,5),  EXXD  (100,5),  EYYD  (100,5), 

EULARY  10 

H  EXYD  (100,5),  SINW  (100,5),  TXX  (100,5),  TYY  (100,5), 

EULARY  11 

I  TXY  (100,5),  SXXSUM(  100,5),  SYYSUM(  100,5),  SXYSUM(100,5), 

EULARY  12 

J  EXEUL1(100,5),  EXEUL2(  100,5),  EXEUL3(  100,5),  NBTYP2(100,5), 

EULARY  13 

K  SMBl  (100,5),  SMB2  (100,5),  EULARYEX 

1 

L  UXl  (100,5),  UYl  (100,5),  UX2  (100,5),  UY2  (100,5), 

EULARYEX  2 

M  XMOMl  (100,5),  YMOMl  (100,5),  XMOM2  (100,5),  YMOM2  (100,5), 
EULARYEX  3 

N  VOLFR  1(100,5),  VOLFR2(  100,5),  TEMPI  (100,5),  TEMP2  (100,5), 

EULARYEX  4 

O  PARDEN(  100,5),  PARNUM(100,5),  DIVl  (100,5),  DIV2  (100,5), 

EULARYEX  5 

P  TMPINT(  100,5),  IFBRN  (100,5)  EULARYEX 

6 

C 

COMMON  /USER/  NTRUSE.  EXTRAS(!00).  LIMUSE 
USER  3 

C  EULARY  15 

DIMENSION  COLBUF(30000)  EULARYEX 

7 

DIMENSION  DTHICK(  100,5)  EULARY 


17 


uuz  uu  u  uu  uuu  uu 


EQUIVALENCE  (X  (1,1),  COLBUF(l)  )  EULARY 

18 

EQUIVALENCE  (EXEUL2(1.1),  DTHICK(l.l)) 

EULARY  19 

THIS  COMMON  BLOCK  CONTAINS  MATERIAL  DATA  FOR  GRAIN  BUR- 
ING 

COMMON  /GRNMSa  LIMMSC,GRNMSC(20) 


DATA  NCYOL/-  l/,NDYNA/0/,TIMEND/14.E-3/,SAMMAX/0. 175/ 

1,  HLSVEL/100./,PHLS/300.E5/,XHLS/0./ 

2,  SPGCN/397. 1 1/,COVOL/0.00 107/,TMPFL/2000./ 

DATA  PCRIT/164.4E5/,UCRIT/n80./ 

GAMM=EXTRAS(2) 


UXEX  =  0. 
UYEX  =  0. 
ROEX  =  0. 
EINEX=  0. 
PEX  =  0. 
UXIEX  =  0. 
UX2EX  =  0. 
UYIEX  =  0. 
UY2EX  =  0. 
SIEIEX  =  0. 
SIE2EX  =  0. 
PDENEX  =  0, 
VFRIEX  =  0. 
VFR2EX  =  0. 
ROIEX  =  0. 
R02EX  =  0. 
IFUOPX  =  0 


SAMACT  =  TMAEB(l) 

FOR  LOW  DENSITIES  THE  VIRIAL  EOS  EQUALS  THE  COVOLUME  EOS 
RHLS  =  L/(SPGCN*TMPFL/PHLS  +  COVOL) 

IF(NCYOL.EQ.NCYCLE)  GOTO  10 
IF(NCYCLE.LE.(NDYNA+1))  XHLS  =  XBAR 
NCYOL  =  NCYCLE 
XHLS  =  XHLS  +  HLSVEL*DLTH 
10  CONTINUE 
C 
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u  u 


MASS  CALCULATIONS 
IF(MASMOM.EQ.l)  THEN 
PEX  =  P(J,IB) 

VFRIEX  =  VOLFRl(J,IB) 

VFR2EX  =  VOLFR2(J,IB) 

ROIEX  =  RH01(J,IB) 

R02EX  =  RH02(J,IB) 

ELSE 

PEX  =  P(J,ID) 

VFRIEX  =  VOLFRl(J,ID) 

VFR2EX  =  VOLFR2(J,ID) 

ROIEX  =  RH01(J,ID) 

R02EX  =  RH02(J,ID) 

ENDIF 

IF(SAMACT.GT.SAMMAX)  GOTO  1 10 
IF(TIME.GT.TIMEND)  GOTO  120 
IF(XHLS.LT.XBAR)  GOTO  130 
C  ROOT  =  AMAX1(0.,  2.*(PHLS/RHLS  -  P(J,IB)/RH01(J,IB))) 
C  IF(ROOT.LE.O.)  GOTO  110 
C  UYIEX  =  SQRT(  ROOT  ) 

C  UY2EX  =  0. 

C  PEX  =  PHLS 
PEX  =  PCRIT 
VFRIEX  =  1.0 
VFR2EX  =  0.0 
ROIEX  =  RHLS 
R02EX  =  1590. 

UYIEX  =  UCRIT 
UY2EX  =  0. 

SIEIEX  =  TMPFL*GRNMSC(2)/GAMM 
SIE2EX=  0. 

GOTO  140 
no  CONTINUE 

C  OPEN(UNIT=99,nLE=STOP.DAT,STATUS=’OLD’) 

C  WRITE(6,111) 

C  111  FORMAT(’SAMMAX  BEREIKT) 

C  CLOSE(UNIT=99) 

C  PRINT  SAMMAX  BEREIKT’ 

GOTO  140 
120  CONTINUE 

C  OPEN(UNIT=99,nLE=STOP.DAT,STATUS=’OLD’) 
WRITE(6,112) 

112  FORMAT(’TIMEND  BEREIKT’) 

C  CLOSE(UNIT=99) 

PRINT  TIMEND  BEREIKT’ 

GOTO  140 
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130  CONTINUE 

C  OPEN(UNlT=99,FILE=STOP.DAT,STATUS=’OLD’) 

C  WRITE(6,113) 

C  113  FORMAT(’XBAR  BEREIKT’) 

C  CL0SE(UN1T=99) 

C  PRINT  XBAR  BEREIKT’ 

GOTO  140 

C  MOMENTUM  CALCULATIONS 
C  IF(MASMOM.EQ.2)  THEN 
C  PEX  =  P(J,ID) 

C  VFRIEX  =  VOLFRl(J,ID) 

C  VFR2EX  =  VOLFR2(JJD) 

C  ROIEX  =  RH01(J,ID) 

C  R02EX  =  RH02(J,ID) 

C  IF(SAMACT.GT.SAMMAX)  GOTO  120 
C  IF(TIME.GT.T1MEND)  GOTO  120 
C  IF(XHLS.LT.XBAR)  GOTO  120 

C  ROOT  =  AMAX1(0.,  2.*(PHLS/RHLS  -  P(J,ID)/RHOl(J,lD))) 

C  IF(ROOT.LE.O.)  GOTO  120 
C  UYIEX  =  SQRT(  ROOT  ) 

C  UY2EX  =  0. 

C  PEX  =  PHLS 

C  ROIEX  =  RHLS 

C  SIEIEX  =  TMPFL*GRNMSC(2)/GAMM 
C  SIE2EX=  0. 

C  120  CONTINUE 

C  ENDIF 

C 

C 

C 

C 

C  TERMINATION  FOR  SUBROUTINE  EXFLOW 
C 

140  CONTINUE 
RETURN 
END 

SUBROUTINE  EXEOTW(AA.BB,CC,DD,TT.EZERO,RHOE,ALPE,ETA,EINTB,PB, 
2L28/J  95 

1  MPN,U,M)  2L28/J  96 

*  <1  *  *  *  *  4c  4c  *  %  «  4c  4c  *  « ifr  4i  «  4t  *  4c  *  *  jfc «  4c  4c  *  4c  4c «  *  %  «  *  *  *  *  :t<  *  *  *  *  * 

c  *  * 

C  ♦  THE  SUBROUTINE  EXEOSQ  ALLOWS  THE  USER  TO  MODEL  A  NON¬ 
STANDARD  * 

C  *  * 

C  *  EQUATION  OF  STATE  * 

C  ♦  * 
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no  no  non 


c  * 
c  * 
c  * 
c  * 
c  * 
c  * 
c  * 


* 

THE  EQUATION  OF  STATE  HAS  THE  FORMAT 
PE  =  (  AA  +  BB*RHOI*EINTE) 

* 


NOTE  THAT  THIS  FORMAT  SHOULD  BE  EQUIVALENT  TO  THE 


C  * 

C  * 
EXTEMP 
C  * 

C 
C 


ENERGY  DEPENDENCE  OF  THE  TEMPERATURE  AS  DEHNED  IN 


COMMON  AJSER/  NTRUSE,  EXTRAS(IOO).  LIMUSE 
USER  3 
C 


COMMON  /NCYVAR/  NCYCLE,  DLTMIN,  NSEXIT,  NUMPLT,  TIMB,  DLTB, 
DLTH  $/NCYVAR  2 

A  ,  DDSTEP,  ITSTEP,  KTSTEP,  NSWRAP,  NUMWAG,  TIME,  DLTE,  DLTFE 
$/NCYVAR  3 

B  ,  SSSTEP,  JTSTEP,  MTSTEP  $/NCYVAR 

4 

C  $/NCYVAR  5 

COMMON  /MSGVAR/  MSG,  IMAX,  NBCMAX,  DTMIN,  ITMIN,  DDMIN, 
IDSTOP,  MSGVAR  2 

1  NZCYC,  IDMAX,  JMAX,  JLIMIT,  MTPRO,  JTMIN,  SSMIN,  LIMSIG, 
MSGVAR  3 

2  NSBLAY.EINTG,  EKING,  EDISG,  XMOMG,  YMOMG,  ZMASG,  NUMSIG, 
MSGVAR  4 

3  EHTRG,  MCHSUB.MASCAL  MSGVAR 

5 


MSGVAR  6 


DATA  NOP/6/ 

VIRIAL  COEFFICIENTS 

DATA  AVIR0/L0E0O/,AVIRl/LO7E-3/,AVIR2/O.5E-6/.AVIR3/0.3E-9/ 

COLD  CONTRIBUTION 
DATA  COLD/0./ 

C 

AVIRI=LlE-3 
AA  =  0. 

BB  =0. 

CC  =  0. 

DD  =  0. 
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TT  =  0. 

RHOGASO  =  0.86728 
GAMMA  1  =  EXTRAS(2) 

C 

IF(ETA.LE.O.)  GOTO  810 

Q  *^t*ilt*ili**il!it:4!*ilcil!ilsilcil!^:ilc^c*>l!^i^i***4cih***itcic^cli,******************* 

c 

C  THE  VIRIALEOS 

C  APROXIMATES  THE  COVOLUME  EOS  VERY  WELL  UP  TO 
C  A  DENSITY  OF  400.00  ,  SAY  , 

C  l./(l.-COVOL*RHO)  =  1.6667 

C  AO  +  Al*RHO  +  A2*RHO**2  +  A3*RHO**3  =  1.4992 

c 

POWO  =  AVIRO 

POWl  =  AVIR1*(RHOGASO*ETA)**1 
POW2  =  AVIR2*(RHOGASO*ETA)**2 
POW3  =  AVIR3*(RHOGASO*ETA)**3 
AA  =  COLD*(ETA  -1.) 

BB  =  (GAMMA  1  -  L)*ETA/(l-POWl) 

CC  =  COLD 

DD  =  (GAMMAl  -  !.)/(  l-POWl)**2 
GOTO  900 
C 

C  ERRORS 
810  CONTINUE 

WRITE(NOP,8I5)  ETA,NCYCLE,IJ,M 
815  FORMAT(lH  .’ERROR  EXEOSQ  ETA=’,E  13.5,415) 

NSEXIT  =  51 
GOTO  900 
C 

820  CONTINUE 

WRITE(NOP,825)  TERM,ETA,NCYCLE,I,J,M 
825  FORMAT(lH  .’ERROR  EXEOSQ  TERM=’,2E  13.5,415) 

NSEXIT  =  52 
GOTO  900 


900  CONTINUE 
RETURN 
END 

SUBROUTINE  TWPHIN(UX1,UX2,UY1,UY2, 
TWPHIN  2 

2  VOLFRl,VOLFR2, 

3  RHOl  ,RH02  ,SIE1  ,SIE2  , 

4  SAMI  ,SAM2  , 


TWPHIN 

TWPHIN 

TWPHIN 
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6  PARDEN, 

TWPHIN 

6 

6  IFBRN 

> 

TWPHIN 

7 

7  VOLMAS, 

TWPHIN 

8 

8 

TWPHIN 

9 

C 

TWPHIN 

10 

TWPHIN 

11 

C  * 

4c 

TWPHIN 

12 

C  * 

SUBROUTINE  TWPHINI  * 

TWPHIN 

13 

C  * 

4c 

TWPHIN 

14 

C  * 

INITIALISES  TWOPHASE  PROPERTIES 

Ik 

TWPHIN 

15 

C  * 

* 

TWPHIN 

16 

C  * 

AT  PURE  START  RUNS  * 

TWPHIN 

17 

C  * 

* 

TWPHIN 

18 

^  4e4e3fc4e3|e3|c]tc3|e3)c:fc4e4e3(c3k4c4e%4(3f(3|e:tc3fc:fc:(e4c3tt4c3ft:fe4c4e3fc:|c]|c4c^3(c]fc3ic3fc3tc:(c%]fe:{e%:4c3fc3(c 

TWPHIN 

19 

C 

TWPHIN 

20 

C 

TWPHIN 

21 

C 

TWPHIN 

22 

TWPHIN 

23 

C  * 

Ik 

TWPHIN 

24 

C  * 

GASPHASE  PROPERTIES  * 

TWPHIN 

25 

C  * 

* 

TWPHIN 

26 

C  * 

RHOl  VOLFRl  UXl  UYl  SAMI  SIEl 

* 

TWPHIN 

27 

C  * 

* 

TWPHIN 

28 

C  * 

* 

TWPHIN 

29 

C  * 

SOLID  PHASE  PROPERTIES  * 

TWPHIN 

30 

C  * 

* 

TWPHIN 

31 

C  * 

RH02  VOLFR2  UX2  UY2  SAM2  PARDEN 

4c 

TWPHIN 

32 

C  * 

IFBRN  SIE2  * 

TWPHIN 

33 

C  * 

4c 

TWPHIN 

34 

^  :|E  I|i  :|i  Ik  ^  I<I  ^  !|I  !<I  !<t  :|<  Ilcifi  ifE  >(< ’ll  Ik  ^  Ik  4c  4cifc  41  41*1(11(1  lit 

TWPHIN 

35 

C 

TWPHIN 

36 

C 

TWPHIN 

37 

^  4c4c4i4tik4i4i4i4>4‘4t4i4t4c4i4tik4iik4i4tik4c4<4<ik4c4c4<4<4t4c4c4t4i4‘4t4c4c4c4t4c4t4c4t4c4<4c4t4c4c 

TWPHIN 

38 

C  * 

4c 

TWPHIN 

39 

C  * 

4c 

TWPHIN 

40 

C  * 

OUTPUT  RHOl  RH02  SAMI  SAM2 

* 
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TWPHIN 
C  * 


C  *  UXl  UX2  *  TWPHIN  42 

C  *  UYl  UY2  *  TWPHIN  43 

C  *  VOLFRI  *  TWPHIN  44 

C  *  VOLFR2  *  TWPHIN  45 

C  *  SIEI  *  TWPHIN  46 

C  ♦  SIE2  *  TWPHIN  47 

C  *  PARDEN  *  TWPHIN  48 

C  *  IFBRN  *  TWPHIN  49 

C  *  *  TWPHIN  50 

********************** 

TWPHIN  51 
C 

COMMON  AJSER/  NTRUSE,  EXTRAS(IOO),  LIMUSE 
USER  3 

C  TWPHIN  52 

C  INLINE  STATEMENT  FUNCTION  TWP¬ 

HIN  53 

MIJBRN(IBTYP,IBN,IBW)  =  IBTYP*IOO  +  IBN*lU  +  IBW 
TWPHIN  54 

C  TWPHIN  55 

C  TWPHIN  56 

C  TWPHIN  57 

C  MUBRN  =  ABC  YES  NO  TWPHIN 

58 

C  ^  SWITCH  FOR  BURNING  IN  THE  PAST  I  0 

TWPHIN  59 

C  TWPHIN  60 


4c  %  «  4c  ]fi  lii  :|i  *  :ft  :)c  4c  *  <1  *  *  :fi  ifi  III  i|>  :|c  ifi  lie  I|f  i|c «  4;  Ik  ^  4c  4c  *  *  *  *  *  *  *  *  *  * 


^  SWITCH  FOR  BURNING  AT  THIS  MOMENT  I  0 

TWPHIN  6 

^  SWITCH  FOR  METHOD  OF  INTERPHASE 
TEMPERATURE  CALCULATION 


C  ^  SWITCH  FOR  B 

TWPHIN  59 
C 

C  ^  SWITCH  FOR  BU 

TWPHIN  61 
C 

C  SWITCH  FOR  M 

TWPHIN  63 

C  TEMPER. 

TWPHIN  64 
C 
C 

C  BURNING  AT  START  ,  IBWAS  =  I 
67 

IFBRN  =  MIJBRN(0,0,0) 

C 

BBTYPE  =  MOD(IFBRN/I00,10) 

70 

IBNOW  =  MOD(IFBRN/I0,I0) 

71 

IBWAS  =  MODaFBRN.IO) 

72 
C 


TWPHIN 

TWPHIN 
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66 

TWPHIN 


TWPHIN 
TWPHIN  69 

TWPHIN 

TWPHIN 


TWPHIN 


TWPHIN 


C  PARTICLE  DENSITY 
lA 

PARDEN  =  1. 

C 

C 

RHOl  =  0.86728 
78 

RH02  =  1570. 

VOLFRl  =  EXTRAS(l) 

80 

VOLFR2  =  1.- VOLFRl 

81 

UXl  =  0. 

UX2  =  0. 

UYl  =  0. 

UY2  =  0. 

SIEl  =  4.1934E5 

86 

SIE2  =  5.8794E5 
87 

SAMI  =  RH01*V0LFR1*V0LMAS 
TWPHIN  88 

SAM2  =  RH02*V0LFR2*V0LMAS 
TWPHIN  89 
C  ALRE  =1.0 


c 

TWPHIN 

90 

c 

TWPHIN 

91 

c 

TWPHIN 

92 

c 

TWPHIN 

93 

RETURN 

TWPHIN 

! 

END 

TWPHIN 

95 

SUBROUTINE  EXSGVG(SGVG,DSDVFR,VFR2,VFR2NW,PDEN2,PDN2NW,RC, 
1  ISWBUR,I,J,M) 

REAL*8  ANGL,D,DX,GSURF,GVOLUM,L.LX,PD,PDM,PI,PX  1  ,PX2,SIDE, 

+  SSUM,TA,TV,VSUM,WEB,WI,WM.WO,X 
INTEGER  n,JJ,NFORMS.NPERF.NCHECK 
LOGICAL  SWITCH 

PARAMETER  (NMAX  =  100  ,  PI  =  3.1415926535897) 

COMMON  /GEOMTR/  WI,WM,WO,L,PD,PDM,NPERF 
COMMON  /FORMS  /  GSURF(NMAX),GVOLUM(NMAX).NFORMS 


TWPHIN 

TWPHIN  75 
TWPHIN  76 
TWPHIN  77 

TWPHIN 

TWPHIN  79 
TWPHIN 

TWPHIN 

TWPHIN  82 
TWPHIN  83 
TWPHIN  84 
TWPHIN  85 
TWPHIN 

TWPHIN 


DATA  SWITCH /.TRUE./ 

DATA  WI  /1.21E-3/,  WO  /1.22E-3/,  L/15.4E-3/,  PDM  /0.64E-3/ 
DATA  NPERF/7/ 


JJ=1 


IF  (SWITCH)  THEN 
SWITCH=.FALSE. 

CALL  TABFRM 
END  IF 

IF  (PDEN2  .GT.  LOE-4)  THEN 
VGRAIN  =  5.213E-7/PDEN2 
ELSE 

VGRAIN=0. 

END  IF 

C  De  numerieke  waarde  is  het  beginvolume  van  de  kruitkorrel 


777  IF  (GVOLUM(JJ+l)  .GT.  VGRAIN  .AND.  JJ  .LT.  NFORMS)  THEN 
JJ=JJ+1 
GOTO  777 
END  IF 

IF  (JJ  .LT.  NFORMS)  THEN 

SGRAIN=DLAGR(VGRAIN,GVOLUM,GSURF,MAX(lJJ-2),MIN(JJ+3, NFORMS)) 
ELSE 

SGRAIN=GSURF(NFORMS)*VGR.\IN/GVOLUM(NFORMS) 

END  IF 

IF  (VGRAIN  .NE.  0.)  THEN 
DSDVFRsSGRAINA^GRAIN 
ELSE 

DSDVFR=0. 

END  IF 

SGVG=VFR2*DSDVFR 

PDN2NW=PDEN2 

VFR2NW=VFR2 

RETURN 

END 

SUBROUTINE  TABFRM 

REAL*8  ANGL,D,DX,GSURF,GVOLUM,L,LX,PD,PDM,PI,PX  1  ,PX2,SIDE, 

+  SSUM,TA.TV,VSUM,WEB,WI,WM,WO^ 

INTEGER  n,JJ,NFORMS,NPERF,NCHECK 
DIMENSION  SIDE(3, 4), ANGL(4,4),NCHECK(4) 

EXTERNAL  GENIS.GENOS 
INTRINSIC  DFLOAT,DMINl,DSQRT 

PARAMETER  (NMAX  =  100  ,  PI  =  3.1415926535897) 

COMMON  /GEOMTR/  WI,WM,WO,UPD,PDM,NPERF 
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COMMON  /FORMS  /  GSURF(NMAX).GVOLUM(NMAX),NFORMS 
DATA  NCHECK/4*-!/ 

PXl  =  0. 

PX2  =  0. 

IF  (PD  .EQ.  0.)  PD  =  PDM 

IF  (NPERF  .EQ.  0)  THEN 
D  =  WI 

WEB  =  DMINl  (D,L) 

ELSE  IF  (NPERF  .EQ.  1)  THEN 
D  =  WI 

WEB  =  DMINl  (0.5*(D-PDM),L) 

ELSE  IF  (NPERF  .EQ.  7)  THEN 
D  =  PDM  +  2.*(WI+WO+PD) 

WEB  =  DMINl  (WI,WO,L) 

SIDE(1,1)  =  WI  +  0.5*(PDM  +  PD) 

SIDE(2,1)  =  SIDE(1,1) 

SIDE(3,1)  =  SIDE(1,1) 

ELSE  IF  (NPERF  .EQ.  19)  THEN 
D  =  PDM  +  4.*PD  +  2.*(WI+WM+WO) 

WEB  =  DMINl  (WI,WM,WO,L) 

SIDE(1,1)  =  WI  +  0.5*(PDM  +  PD) 

SIDE(2,1)  =  SIDE(Ll) 

SIDE(3,1)  =  SIDE(1,1) 

SIDE(1,2)  =  0.5*DSQRT(3.*(WM+PD)**2  +  (WI+0.5*(PDM+PD))**2) 

SIDE(2,2)  =  SIDE(1,2) 

SIDE(3,2)  =  SIDE(1,1) 

SIDE(1,3)  =  0.5*(WI+WM)  +  0.25*PDM  +  0.75*PD 
SIDE(2,3)  =  SIDE(1,2) 

SIDE(3,3)  =  WM  +  PD 
SIDE(1,4)  =  SIDE(1,3) 

SIDE(2,4)  =  2.*SIDE(1,3) 

SIDE(3.4)  =  SIDE(L3)*DSQRT(3.0D0) 

END  IF 

Vul  arrays  met  brandend  oppervlak  en  resterend  volume 
i  =  50  (=nmax/2)  bij  doorbranden  van  web  (brand  aan  twee  kanten) 

DO  1  II  =  LNMAX 

Bereken  ingebrande  geometrie  variabelen 
X  =  DFLO  AT(II- 1  )*  WEB/NMAX 
LX  =  L  -  2.*X 
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DX  =  D  -  2  *X 
SSUM  =0. 

VSUM  =0. 

IF  (NPERF  .GE.  1)  PXl  =  PDM  +  2  *X 
IF  (NPERF  .GE.  7)  PX2  =  PD  +  2.*X 


Als  web  nog  niet  doorgebrand  dan  directe  berekening 
IF  (2.*X  .LE.  WEB)  THEN 

TA  =  0.25*PI*(DX**2  -  PX1**2  -  (NPERF-1)*PX2**2) 

SSUM  =  2.*TA  +  LX*PI*(DX  +  PXI  +  (NPERF- 1)*PX2) 

VSUM  =TA*LX 

Web  doorgebrand  dan  verder  rekenen  (alleen  7  en  19  gats) 

ELSE  IF  (NPERF  .EQ.  7)  THEN 

CALL  GENIS  (SIDE(1,1),ANGL(1,1),PX1,PX2,LX,NCHECK(1),TA,TV) 
SSUM  =SSUM  +6.*TA 
VSUM  =VSUM  +6.*TV 

CALL  GENOS  (SIDE(1,1),ANGL(1,1).PX2,LX,DX/2.,NCHECK(2),TA,TV) 
SSUM  =SSUM  +6.*TA 
VSUM  =VSUM  +6.*TV 

ELSE  IF  (NPERF  .EQ.  19)  THEN 

CALL  GENIS  (SIDE(1.1),ANGL(1,I),PX1,PX2,LX,NCHECK(1),TA,TV) 
SSUM  =SSUM  +6.*TA 
VSUM  =VSUM  +6.*TV 

CALL  GENIS  (SIDE(1,2),ANGL(1,2).PX2,PX2,LX,NCHECK(2),TA,TV) 
SSUM  =SSUM  +6*TA 
VSUM  =VSUM  +6.*TV 

CALL  GENIS  (SIDE(1,3),ANGL(1,3),PX2,PX2,LX,NCHECK(3),TA,TV) 
SSUM  =  SSUM  +  12.*TA 
VSUM  =  VSUM  +  12.*TV 

CALL  GENOS  (SIDE(1,4),ANGL(1.4),PX2,LX.DX/2.,NCHECK(4),TA,TV) 
SSUM  =  SSUM  +  12.*TA 
VSUM  =  VSUM  +  12.*TV 

END  IF 

GSURF(n)  =SSUM 
GVOLUMai)  =  VSUM 
IF  (VSUM  .EQ.  0.)  GOTO  2 

1  CONTINUE 

2  NFORMS  =  11-1 

RETURN 

END 
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*  * 
* 

* 

* 

* 

* 

* 

* 

* 

*  * 


i^ilnif^iifltlf^iflti^ik********************** 

* 

Function  :  DLAGR  * 

Datum  Vl.O  :  04  januari  1989  * 

* 

Ondenverp  :  Interpolatie  * 

Doel  :  Interpoleren  met  methode  van  Lagrange  * 

Auteur  :  J.C.  Post  * 

* 


DOUBLE  PRECISION  FUNCTION  DLAGR(X^TAB.YTAB,ELMNTI,ELMNT2) 

REAL*8  X,XTAB.YTAB,PROD 
INTEGER  n,JJ,ELMNTl,ELMNT2 
DIMENSION  XTAB(*),YTAB(*) 

DLAGR  =0.  D-KX) 

DO  2  II  =  ELMNT1,ELMNT2 
PROD  =  1.  D+00 
DO  1  JJ  =  ELMNT1,ELMNT2 

IF  (JJ  .NE.  H)  PROD  =  PROD*(X-XTAB(JJ))/(XTAB(II)-XTAB(JJ)) 

1  CONTINUE 

DLAGR  =  DLAGR  +  YTAB(n)*PROD 

2  CONTINUE 

RETURN 

END 

i^^^citiilit****illL******i^***^^1^*1f1H^****^^*******i^*********ii^**1^******^^*********i^** 


SUBROUTINE  GENIS  (S.A,Rl,R2,LX,IFLAG,SURF,VOL) 


*  subroutine  *genis*  :  calculate  surface  area  and  volume  for  a 

*  general  inner  sliver  of  a  burning  grain  with  length  =  LX, 

*  mid  perf  diam  =  R1  and  other  perf  diam  =  R2. 

REAL*8  A,B2,B3,E,LX,PI,RAD,R  1  ,R2,S,SURF,TAU  1 2,TAU  1 3,TAU2 1 , 

+  TAU23,TAU31,TAU32,VOL 
INTEGER  n,JJ,IFLAG 

DIMENSION  S(3),A(4) 

DATA  PI  /  3.1415926535897  / 

VOL  =  0. 

SURF  =0. 

*  initial  pass  :  iflag  was  set  negative  by  calling  routine. 
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*  store  angles  al,a2,a3  and  area  of  triangle 

*  with  sides  s(l),s(2),s(3)  into  a(l) . a(4). 

IF  (IFLAG  .LT.  0)  THEN 
IFLAG  =0 

A(l)  =  ACOS((S(2)**2+S(3)**2-S(l)**2)/(2  *S(2)*S(3))) 
A(2)  =  ACOS((S(l)**2+S(3)**2-S(2)**2)/(2  *S(1)*S(3))) 
A(3)  =PI-A(1)-A(2) 

A(4)  =  0.5*S(1)*S(3)*SIN(A(2)) 

*  check  for  error  condition  :  find  if  triangle  acceptable. 

JJ  =0 

DO  15  n  =  1,3 

IF  (Aai)  -LT.  0.25*PI)  JJ  =  JJ  +  1 
15  CONTINUE 

IF  (JJ  .GT.  1)  THEN 

PRINT  Error  at  subroutine  *GENIS*  (ANGLE  <  PI/4)  ’ 
STOP 
END  IF 

END  IF 

*  succeeding  passes  until  burnout :  find  auxiliary  angles. 

B2  =  0.5*(S(2)  +  0.25*(Ri**2  -  R2**2)/S(2)) 

B3  =  0.5*(S(3)  +  0.25*(R1**2  -  R2**2)/S(3)) 

TAU12  =  ACOS(DMINI(1.0D0,B3/(0.5*RI))) 

TAU21  =  ACOS(DMINl(1.0D0,(S(3)-B3)/(0.5*R2))) 

TAU13  =  ACOS(DMINl(1.0D0,B2/(0.5*Rl))) 

TAU31  =  ACOS(DMINl(1.0D0,(S(2)-B2)/(0.5*R2))) 

TAU23  =  ACOS(DMINl(1.0D0,S(l)/R2)) 

TAU32  =TAU23 


*  and  branch  if  sliver  fails  burnout  criteria. 

IF  (TAU12+TAU13+TAU23+TAU21+TAU31+TAU32  .LT.  PI  .AND. 
#  LX  .GT.  0.)  GOTO  25 

*  sliver  just  burned  out :  set  flag  to  bypass  area  and  volume  calculati 

IFLAG  =  1 
GOTO  30 


*  sliver  not  burned  out :  determine  end  area  ,  volume  and  surface  area. 

25  E  =  A(4)  -  0.25*R2*(  S(1)*SIN(TAU23)  +  S(3)*SIN(TAU12)  + 
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#  S(2)*S1N(TAU31)  )  - 

#  (A(l)-TAU13-TAU12)/8.*Rl**2  - 

#  (A(2)-TAU21-TAU23  +  A(3)-TAU31-TAU32)/8.*R2**2 
VOL  =  E*LX 

SURF  =  2  *E  +  0.5*LX*R1*(PI/3.-TAU12-TAU13)  + 

#  0.5*LX*R2*(2./3*PI-TAU21-TAU23-TAU31-TAU32) 


*  and  return  to  main 

RETURN 

*  sliver  is  burned  out  :  return  with  zero  volume  and  surface  area. 

30  VOL  =  0. 

SURF  =0. 

RETURN 

END 

SUBROUTINE  GENOS(S,A,R2,LX,RAD,IFLAG.SURF,VOL) 


*  subroutine  *genos*  ;  calculate  surface  area  and  volume  for  a 

*  general  outer  sliver  of  a  burning  grain 

*  with  length  =  LX  ,  radius  =RAD  and 

*  perf  diam  =  R2.  ,  radius  =RAD  and 

REAL*8  A,E,LX,RAD,R2,S,SIG,SURF,TAU1,TAU2,TAU3,TAU4,V0L 
INTEGER  IFLAG 

DIMENSION  S(3),A(4) 

IF  (IFLAG)  10,20,30 

*  initial  pass  :  iflag  was  set  negative  by  calling  routine, 

*  store  angles  al,a2,a3  and  area  of  triangle 

*  with  sides  s(l),s(2),s(3)  into  a(I) . a(4). 

10  A(l)  =  ACOS((S(2)**2+S(3)**2-S(l)**2)/(2.*S(2)*S(3))) 

A(2)  =  ACOS((S(l)**2+S(3)**2-S(2)**2)/(2.*S(l)*S(3))) 

A(3)  =  ACOS((S(l)**2+S(2)**2-S(3)**2)/(2.*S(l)*S(2))) 

A(4)  =  0.5*S(1)*S(3)*SIN(A(2)) 

*  set  flag  to  zero  to  bypass  initialization  hereafter. 
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IFLAG  =  0 


succeeding  passes  until  bumout :  first  determine  auxiliary  angles. 


20  TAUl  =  ACOS(DMINl(1.0D0,(S(2)**2+RAD**2-0.25*R2**2)/ 

#  (2.*S(2)*RAD))) 

TAU2  =  ACOS(DMINl(1.0D0,(S(3)**2+RAD**2-0.25*R2**2)/ 

#  (2.*S(3)*RAD))) 

TAU3  =  ACOS(DMAXl(-1.0D0,(S(2)**2-RAD**2+0.25*R2**2)/ 

#  (S(2)*R2))) 

TAU4  =  ACOS(DMAXU-1.0D0,(S(3)**2-RAD**2-K).25*R2**2)/ 

#  (S(3)*R2))) 

SIG  =  ACOS(DMIN1(1.0DO,S(1)/R2)) 


then  check  error  conditions 

IF  (TAU3  .LT.  A(3)  .OR.  TAU4  .LT.  A(2))  THEN 
STOP  ’  ERROR  AT  SUBROUTINE  *GENOS*  ’ 
END  IF 


if  ok,  branch  if  sliver  fails  bumout  criteria. 

IF  (TAU1+TAU2  .LT.  A(l)  .AND.  LX  .GT.  0.)  GOTO  25 

sliver  just  burned  out :  set  flag  to  bypass  area  and  volume  calculati 

IFLAG  =  1 
GOTO  30 


sliver  not  burned  out :  determine  end  area  ,  volume  and  surface  area. 

25  E  =  0.5*RAD*(S(2)*SIN(TAU1)+RAD*(A(1)-TAU1-TAU2)  + 

#  S(3)*SIN(TAU2))  -  A(4)  -  0.25*R2*(S(1)*SIN(SIG)  + 

#  0.5*R2*(TAU3+TAU4-2.*SIG-A(2)-A(3))) 

VOL  =  E*LX 

SURF  =  2.*E  +  LX*(RAD*(A(1)-TAU1-TAU2)  + 

#  0.5*R2*(TAU3+TAU4-2.*SIG  -  A(2)  -  A(3))) 


and  return  to  main 
RETURN 


sliver  is  burned  out :  return  with  zero  volume  and  surface  area. 

30  VOL  =  0. 

SURF  =0. 


RETURN 


END 


SUBROUTINE  EXEDIT(NCOL.NBC) 

4c  *  :|c «  4c  *:(<  :|ci|c  *4:  >1:  *  !|<  ****  If:  i<<  lie  *  :fE  4ciiE  ilciii  ik  ifi  %******  ******  *****:!<**  i|c  **  i)c  * 

*  * 

*  THIS  SUBROUTINE  ALLOWS  THE  USER  ADDITIONAL  EDITING 


*  12-4-91  M.  JURGENS 

*  * 

*****************************c(,*:|i4,:|c******************  *************** 


LOGICAL  SSWITCH 


THIS  COMMON  BLOCK  CONTAINS  MATERIAL  DATA  FOR  GRAIN  BUR- 


iNllNLj 


COMMON  /GRNMSC/  LIMMSC,GRNMSC(20) 


COMMON  AJSER/  NTRUSE,  EXTRAS(IOO),  LIMUSE 
USER  3 

C  USER  4 

COMMON  /CUTOFF/  NTRCUT,  VELCUT,  VELMAX,  RHOCUT,  SSLLIM, 
FBLEND,  CUTOFF  2 

1  SSLLSQ,  VFRCUT  CUTOFF  3 

C  CUTOFF  4 

C  BDARY  8 

C  THIS  DIMENSION  STATEMENT  ASSUMES  L1MROW(2)=100  AND  NBCOLS(2- 
)=5  EULARY  2 

COMMON  /COLMEM/  LEMCOB,  EULARY 

3 

AX  (100,5),  Y  (100,5),  THICK  (100,5),  UX  (100,5),  EULARY 

4 


B  UY  (100,5),  XMOM  (100,5),  YMOM  (100,5),  VOLMAS(  100,5), 

EULARY  5 

C  VOID  (100,5),  RHOl  (100,5),  RH02  (100,5),  SAMI  (100,5), 

EULARY  6 

D  SAM2  (100,5),  SIEl  (100,5),  SIE2  (100,5),  P  (100,5),  EULARY 

7 

E  Q  (100,5),  SSPD  (100,5),  DIV  (100,5),  ERGl  (100,5),  EULA¬ 

RY  8 

F  ERG2  (100,5),  NBTYPE(  100,5).  NGATES(  100,5),  FCOVV  (100,5), 

EULARY  9 

G  FCOVSK  100,5),  FCOVSJ(100,5),  EXXD  (100,5),  EYYD  (100,5), 

EULARY  10 

H  EXYD  (100,5),  SINW  (100,5),  TXX  (100,5),  TYY  (100,5), 
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EULARY  1 1 

I  TXY  (100,5),  SXXSUM(100,5),  SYYSUM(  100,5),  SXYSUM(I00,5), 

EULARY  12 

J  EXEUL1(100,5),  EXEUL2(100,5),  EXEUL3(  100,5),  NBTYP2(  100,5), 

EULARY  13 

K  SMBl  (100,5),  SMB2  (100,5),  EULARYEX 

1 

L  UXl  (100,5),  UYl  (100,5),  UX2  (100,5),  UY2  (100,5), 

EULARYEX  2 

M  XMOMl  (100,5),  YMOMl  (100,5),  XMOM2  (100,5),  YMOM2  (100,5), 
EULARYEX  3 

N  VOLFR1(100.5),  VOLFR2(  100,5),  TEMPI  (100,5),  TEMP2  (100,5), 

EULARYEX  4 

O  PARDEN(100,5),  PARNUM(  100,5),  DIVl  (100,5),  DIV2  (100,5), 

EULARYEX  5 

P  TMPINT(  100,5),  IFBRN  (100,5)  EULARYEX 

6 

C  EULARY  20 

COMMON  /NCYVAR/  NCYCLE,  DLTMIN,  NSEXIT,  NUMPLT,  TIMB,  DLTB, 
DLTH  $/NCYVAR  2 

A  ,  DDSTEP,  ITSTEP,  KTSTEP,  NSWRAP,  NUMWAG,  TIME,  DLTE,  DLTFE 
$/NCYVAR  3 

B  ,  SSSTEP,  JTSTEP,  MTSTEP  $/NCYVAR 

4 

C  $/NCYVAR  5 

COMMON  /COLVAR/  NCOLC,  EINTC(5),  XMOMC(5),  lA,  ID, 

COLVAR  2 

1  NCOLA,  NCOLD,  EKINC(5),  YMOMC(5),  IB,  IE, 

COLVAR  3 

2  NCOLB,  NCOLE,  EDISC(5),  ZMASC(5),  EHTRH(5).  IC 

COLVAR  4 

C  COLVAR  5 

COMMON  /MSGVAR/  MSG,  IMAX,  NBCMAX,  DTMIN,  ITMIN,  DDMIN, 
IDSTOP,  MSGVAR  2 

1  NZCYC,  IDMAX,  JMAX,  JLIMIT,  MTPRO,  JTMIN,  SSMIN,  LIMSIG, 

MSGVAR  3 

2  NSBLAY,EINTG,  EKJNG,  EDISG,  XMOMG,  YMOMG,  ZMASG,  NUMSIG, 
MSGVAR  4 

3  EHTRG,  MCHSUB,MASCAL  MSGVAR 

5 

C  MSGVAR  6 

COMMON  /RUNVAR/  NCYLIM,  NCYREF,  NSYM,  GRAVX,  LIMPAG,  DATE(5), 
RUNVAR  2 

1  HEAD(18),  TIMLIM,  ENFRAC,  FRAT,  GRAVY,  LINEPP,  HOUR(2), 

RUNVAR  3 

2  NCYBEG,  DLTWRP,  STATIC,  KENT,  XBDAMP,  YBDAMP,  EDAMP, 
RUNVAR  4 

3  METHDT,  DLTMAX,  MAXSCM,  MAXLCM,  MAXDSK,  PHEAD(8), 


68 


NDELK,  RUNVAR  5 

4  NILSUP,  IFEVPL,  IFEVPE,  NSYMDF  RUNVAR 

6 

C  RUNVAR  7 

C 


DATA  SSWITCH/.TRUE./ 

DATA  NPRINT/100/ 

IF  (SSWITCH)  THEN 
PPMAX=0. 

IVAR=0 

JVAR=0 

MVAR=0 

TVAR=0. 

DO  nil  JJZ=1,100 
DO  1112IIZ=1,5 
EXEUL1(JJZ,IIZ)=0. 

EXEUL2(JJZ,IIZ)=0 

EXEUL3(JJZ,IIZ)=0. 

1112  CONTINUE 
nil  CONTINUE 

SSWITCH=.FALSE. 

ENDIF 

DO  no  JZ  =  2,JMAX 

IF(SAM1(JZ,NBC).LE.0.  .AND.  SAM2(JZ,NBC).LE.O.)  GOTO  110 
IF  (P(JZ,NBC)  .GT.  PMAXX)  THEN 
PMAXX=P(JZ,NBC) 

IVAR=NCOL 

JVAR=JZ 

MVAR=NCYCLE 

TVAR=TIME 

ENDIF 

IF  (P(JZ,NBC)  .GT.  EXEUL1(JZ.NBC))  THEN 
EXEULl  (JZ,NBC)=P(JZ,NBC) 

EXEUL2(JZ,NBC)=NCYCLE 

EXEUL3(JZ,NBC)=TIME 

ENDIF 

no  CONTINUE 

IF  (MOD(NCYCLE,NPRINT).EQ.O)  THEN 
WRITE(6,9388)TVAR,MVAR,IVAR,JVAR.PMAXX 
9388  FORMAT(lH  ,’**INFO  EXEDIT,  MAXIMUM  PRESSURE\E13.5,3I6,3X.E13.5) 
ENDIF 
C 

C  EXTRAS(21)  AND  (22)  ARE  DEFINED  IN  EXBDFO  **************** 
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c 

EXEUL1(1,NBC)=EXTRAS(21) 

EXEUL2(  1  .NBC)=EXTRAS(22) 

EXRIG1=EXTRAS(23) 

EXRIG2=EXTRAS(24) 

EXRIG3=EXTRAS(25) 

C  WRITE(NOP,740)  EXRIG  1,EXRIG2,EXRIG3 
C740  F0RMAT(/,8X,’ 

C  1  EXRIGl(23)’,E14.5y8X.’EXRIG2(24)=’,E14.5y8X,’ 

C  2  EXRIG3(25)=\E14.5/) 

C 

C 

RETURN 

END 

Q  **tt***^t*ii^c*^i*^iif**if*if^i*^L^i*ici^:tticif^c*ik**************************************** 


SUBROUTINE  EXBDFO 
C 

C  THIS  COMMON  BLOCK  ASSUMES  LIMBD  =  500 
BDARY  2 

COMMON  /BDARY/  NTRBDY,  XENDBD(500),  XDOTBD(500).  XFORBD(500), 
BDARY  3 

1  NTYPBD(500),  YENDBD(500),  YDOTBD(500),  YFORBD(500), 

BDARY  4 

2  PMASBD(500),  DTDXBD(500),  DTDYBD(500),  AREABD(500), 

BDARY  5 

3  WORKVI,  XIMPVI,  YIMPVI,  NUMBD,  LIMBD, 

BDARY  6 

4  DTDXBS(500),  DTDYBS(500)  BDARY 

7 

C 

COMMON  /NCYVAR/  NCYCLE,  DLTMIN,  NSEXIT,  NUMPLT,  TIMB,  DLTB, 
DLTH  $/NCYVAR  2 

A  ,  DDSTEP,  ITSTEP,  KTSTEP,  NSWRAP,  NUMWAG,  TIME,  DLTE,  DLTFE 
$/NCYVAR  3 

B  ,  SSSTEP,  JTSTEP,  MTSTEP  $/NCYVAR 

4 

C  $/NCYVAR  5 

COMMON  /MSGVAR/  MSG,  IMAX,  NBCMAX,  DTMIN,  ITMIN,  DDMIN, 
IDSTOP,  MSGVAR  2 

1  NZCYC,  IDMAX,  JMAX,  JLIMIT,  MTPRO,  JTMIN,  SSMIN,  LIMSIG, 
MSGVAR  3 

2  NSBLAY,EINTG,  EKING,  EDISG,  XMOMG,  YMOMG,  ZMASG,  NUMSIG, 
MSGVAR  4 

3  EHTRG,  MCHSUB,MASCAL  MSGVAR 

5 

C  MSGVAR  6 
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n  n  c  con 


THIS  COMMON  BLOCK  ASSUMES  LIMUSE  =  100 
SER  2 

COMMON  AJSER/  NTRUSE,  EXTRAS(IOO),  LIMUSE 
SER  3 

USER  4 


DIMENSION  LNTPRJ(4) 

C 

DATA  LNTPRJ/1,2,3,4/,  NSWOFF/0/,  IPl/l/,  ARPRO/4.53646E-3/ 

DATA  NUMPRJ/4/,  POFF/5.E6/,  PATM/1.E5/.  XPO/0.588/ 

DATA  XD1^D2,XD3^D4/.002,.012,.024,.040/ 

DATA  PFl/0./ 

DATA  NOP/6/,  NCALL/0/ 

C 

PF2=EXTRAS(3) 

PF4=EXTRAS(4) 

IF  (NCALL.EQ.0)  THEN 
FACO  =  (PF1-P0FF)/XD1 
FACl  =  (PF2-PF1)/(XD2-XD1) 

FAC2  =  (PF4-PF2)/(XD4-XD3) 

ENDIF 

NCALL  =  NCALL  +  1 
FXTOFF  =  ARPRO  *  POFF 
C 

FXTOT  =  0. 

DO  10  JJ  =  l.NUMPRJ 
II  =  LNTPRJ(JJ) 

FXTOT  =  FXTOT  +  XFORBD(II) 

10  CONTINUE 
C 

IF  (NSWOFF  .EQ.  0)  THEN 
IF  (FXTOT  .GT.  FXTOFF)  THEN 
C  THE  PROJECTILE  WILL  START  TO  MOVE 

NSWOFF  =  1 

WRITE(NOP,610)  TIME,FXTOT,FXTOFF 
610  FORMAT(/,8X,’EXDBFO:  THE  PROJECTILE  STARTS  TO  MOVE  AT  T 
1  E14.5y8X,’  FXTOT=’,E14.5,’  FXOFF=’,E14.5/) 

ENDIF 

ENDIF 

C 

IF  (NSWOFF  .EQ.  0)  THEN 
C  NSWOFF  =  0  :  SET  FORCES  TO  ZERO. 

DO  20  JJ  =  LNUMPRJ 
II  =  LNTPRJ(JJ) 

XFORBD(n)  =  0. 

20  CONTINUE 
ELSE 


XDIS  =  XENDBD(LNTPRJ(1))  -  XPO 
IF  (XDIS  .LT.  XDl)  THEN 
PFRIC  =  POFF  +  XDIS*FAC0 
ELSE  IF  (XDIS  .LT.  XD2)  THEN 
PFRIC  =  PFl  +  (XDIS-XD1)*FAC1 
ELSE  IF  (XDIS  .LT.  XD3)  THEN 
PFRIC  =  PF2 

ELSE  IF  (XDIS  .LT.  XD4)  THEN 
PFRIC  =  PF2  +  (XDIS-XD3)*FAC2 
ELSE 

PFRIC  =  PF4 
ENDIF 
C 

FXFRIC  =  ARPRO*(PFRIC  +  PATM) 

IF  (XDOTBD(IPl)  .GT.  0.) 

*  XFORBDaPl)  =  XFORBDaPl)  -  FXFRIC 
ENDIF 
C 

EXTRAS(21)  =  XDIS 
EXTRAS(22)  =  PFRIC 
C 

EXTRAS(23)  =  FXTOT 


EXTRAS(24)  =  FXTOT-FXFRIC 
EXTRAS(25)  =  FXTOT/ARPRO 

C  WRITE(NOP,730)  EXTRAS(23),EXTRAS(24),EXTRAS(25) 

C  730  FORMAT(/,8X,’blablabla 

C  1  EXTRAS(23)’,E14.5,/8X,’EXTRAS(24)=’,E14.5y8X,’ 

C  2  EXTRAS(25)=’,E14.5/) 

C 

c 

RETURN 

END 

C 

^il‘iii‘ic;liifiilt**’li*i‘**>l‘*’l<***il‘i<****^**’l‘i‘**ilii‘*!l‘!l‘**i‘**********!li*>li***** 

c 

SUBROUTINE  EXBUR(BURN,DBDRH0,DBDVFR,RH0l,RH02, 

1  VFR1,VFR2,UXI,UYI,UX2.UY2. 

2  SIE1,SIE2,PDEN2,PB, 

3  I.J,M) 

C 

DATA  CONST/.Ol/  BETA/8.E-10/ 


BURN=BETA*PB+CONST 

DBDRHO=0. 

DBDVFR=0. 

RETURN 

END 
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